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CN i ABSTRACT 

We present results from a comprehensive number of relativistic, time-dependent, ax- 
isymmetric simulations of the runaway instability of non-constant angular momentum 
I thick discs around black holes. This second paper in the series extends earlier results 

where only constant angular momentum discs were considered. All relevant aspects of 
the theory of stationary thick discs around rotating black holes, necessary to build the 
equilibrium initial data used in our simulations, are presented in great detail. The an- 
gular momentum of the evolved discs is assumed to increase outwards with the radial 
I distance according to a power law, I — Kr" , where K > corresponds to prograde 

. discs (with respect to the black hole rotation) and K < to retrograde discs. The 

main simplifying assumptions of our approach are not to include magnetic fields and 
self-gravity in the discs (test-fluid approximation). Furthermore, the dynamics of the 
spacetime is accounted for by computing the transfer of mass and angular momentum 
O ' from the disc to the black hole through the event horizon. In this approximation the 

evolution of the central black hole, which initially is non-rotating, is assumed to follow 
^ • a sequence of Kerr black holes of increasing mass and spin. All discs we build slightly 

^ ' overflow the potential barrier at the cusp, departing from equilibrium, so that accre- 

tion is possible. In agreement with previous results based on stationary models we find 
that by allowing the mass and the spin of the black hole to grow, constant angular 
rS I momentum discs rapidly become unstable on a dynamical timescale (a few orbital pe- 

^ '. riods). The comparison with the results of our first paper shows that the effect of the 

angular momentum transfer from the tori to the black hole is to make constant angular 
momentum discs less unstable, increasing the timescale for the runaway instability to 
grow. However, we find that non-constant angular momentum discs are dramatically 
stabilized for very small values of the angular momentum slope a, much smaller than 
the Keplerian value a = 1/2. Our fully relativistic and time-dependent simulations 
confirm, thus, the predictions of stationary studies concerning the stabilizing effect 
of non-constant angular momentum distributions. For the various disc-to-hole mass 
ratios considered we systematically find that the critical values of a below which the 
runaway instabilty can exist are slightly smaller than those reported previously in the 
literature based on stationary studies. 

Key words: accretion, accretion discs - black hole physics - hydrodynamics - insta- 
bilities - gamma rays: bursts. 
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1 INTRODUCTION 

In a recent paper jFont fc Daignd l2002lJ) (Paper I here- 
after) we presented the first time-dependent, hydrodynami- 
cal simulations in general relativity of the so-called runaway 
instability of geometrically thick discs (or tori) around a 
Schwarzsc hild black hole. The origin of this instability, dis- 
covered bv lAbramowicz et al.l il983l) . is a dynamical process 



by which the cusp of a disc which initially is filling its Roche 
lobe penetrates inside the disc due to mass transfer from the 
disc to the accreting black hole. This process leads to the 
complete destruction of the disc on a dynamical timescale. 
lAbramowicz et all (^8^ found that the runaway instability 
occurs for a large range of parameters such as the disc-to- 
hole mass ratio and the location of the inner radius of the 
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disc. More detailed studies followed (see, e.g. Paper I for 
an up-to-date list), most of which were based on stationary 
models in which a fraction of the mass and angular momen- 
tum of the initial disc is instantaneously transferred to the 
black hole. The new gravitational field allows to compute the 
new position of the cusp, which controls the occurrence of 
the runaway instability. The self-gravity of the tori favours 
the instability, as shown from bo th, studies based on a 
aseudo-potential for the b lack hole (|Khanna fc Chakrabartil 



Il992l : iMasuda et alJll998h and from fully relativistic (sta- 
tionary) calculations ll^shida et aDliggGl) . However, there 
are some parameters in the models which have st abilizing 
effects. One of this is the rotation of the black hole JWilsoni 
ll984l : lAbramowicz et al.lll998ll . but the most important one 
seems to be the radial distribution of the specific angular 
momentum, increasing with the radial distance as a power 
law (Daigne & Mochkovitch 1997; Abramowicz ct al. 1998). 
Despite the abundant literature on the subject, it is com- 
monly accepted that the very existence of the instability 
remains still uncertain. 

In Paper I we initiated a program aimed at investigat- 
ing the runaway instability through time-dependent hydro- 
dynamical simulations of increasing complexity in general 
relativity. The starting point was simple. In the simulations 
reported in Paper I the distribution of the angular momen- 
tum in the disc was assumed to be constant. We also ne- 
glected the self-gravity of the disc as well as the presence 
of magnetic fields. Our model, however, accounted for the 
dynamics of the gravitational field of the black hole plus 
torus system in an approximate yet satisfactory manner. 
The evolution of the central black hole was assumed to fol- 
low a sequence of Schwarzschild black holes of increasing 
mass, whose increase rate was controlled by computing the 
mass flux from the disc to the black hole through the inner- 
most point of the numerical grid. We notice that the ability 
to take into account such dynamics is of paramount impor- 
tance to investigate the runaway instability at all. 

Our study showed that by allowing the mass of the black 
hole to grow the runaway instability appears on a dynam- 
ical timescale, in agreement with previous estimates from 
stationary models. Since our study was mainly motivated 
by understanding the hydrodynamical evolution of the cen- 
tral engine of gamma-ray bursts (GRBs hereafter), we only 
considered the case of stellar mass black holes, leaving aside 
thick discs around supermassive black holes which can also 
be present in active galactic nuclei. For a black hole of 2.5Mq 
and disc-to-hole mass ratios in the range 0.05 to 1, we found 
in Paper I that the timescale of the instability never exceeds 
1 s for a large range of mass fluxes and it is typically a few 10 
ms. We note in passing that the runaway instability seems 
to be a robust feature of non-self-gravitating, constant an- 
gular momentum discs irrespective of the way accretion is 
induced, aSiZanotti ct al. (2003) have recently reported. 

Whatever the progenitor system which collapses to 
form a stellar mass black hole surrounded by a thick disc 
is, it is very unlikely that the initial distribution of the 
angular momentum in the disc is constant. As numeri- 
cal simulations show this is not the case in discs formed 
after the coalescence of a close compact binary (consist- 
ing either of two neutron stars or of a black hole and 



a neutron star) (|Dg^ie^^^J_^2ijBiUfi££t-£^iJ 

iKluzniak fc Led Il998l: iRuffert fc Jankal Il999l: iLe 



200C; 



IShibata et alj|2 003[). in the coUapsa r model of failed super- 
nova llMacFadven fc Wooslevlll999l) . or, at different scales, 
in the collapse of a rotating supermassive star leading to a 
rotating superm assive black hole tShibata fc Shapiro 200*3 : 
IShapiro fc Shib ata 2002 :). For instance, in the Newtonian 
SPH simulations of yDavies et alj il994r) the final configu- 
ration consists of a core of mass 2.44M0 in uniform rota- 
tion surrounded by a differentially rotating thick disc of 
mass O.36M0. The distribution of the angular momentum 
in the disc is close to a power law of index 0.2. Similarly, 
in the mos t recent thr ee-dimensional SPH coalescence sim- 
ulations of iLeel i200(]f) the corresponding angular momen- 
tum distribution follows a power law of index 0.45, quite 
close to Keplerian. We note in particular that recent rel- 
ativisti c simulations of un equal mass neutron star coales- 
cence jShibata et al.' "200^ yield a black hole surrounded 
by a thick disc, followi ng tid al disruption of the less mas- 
sive star. lshibata et al.l i2003i) found that for a neutron star 
rest-mass ratio ~ 0.85 the final mass of the disc may be 
several percents of the total mass of the system if the stars 
are not very compact. Furthermore, the angular momentum 
distribution in the disc that forms is not constant (M. Shi- 
bata, private communication). In addition to the evidence 
provided by the results of numerical simulations, as the dy- 
namical timescale on which the runaway instability occurs is 
much shorter than the viscous timescale necessary to achieve 
angular momentum transport, it seems much more realistic 
to assume a non-constant angular momentum distribution 
in the disc, increasing outwards from the cusp to larger radii. 

Therefore, in the present paper we extend the study 
initiated in Paper I to such most interesting case of non- 
constant angular momentum discs. The main motivation 
of our work is to check through time-dependent simula- 
tions in general relativity whether such distributions have 
indeed the stabilizing effect previou sly found in non self- 
gravitating stationary models (Daigne fc Mochkovitchll997l : 
lAbramowicz et al]ll998l : ILu et aL.2000ih . We note that pre- 
liminary results of t his investigation were a lready presented 
m a previous letter llFont fc Daigndl2002a^ (Paper II here- 
after). In the current paper we extend those results in two 
main directions: firstly, we present in great detail all rele- 
vant aspects of the theory of stationary, non-constant an- 
gular momentum thick discs around rotating black holes, 
which is necessary to build the equilibrium initial data we 
use in our simulations. Secondly, we report on the results of 
a comprehensive number of simulations which allows us to 
explore the dependence of the absence/presence of the run- 
away instability on various parameters of the model such as 
the disc-to-hole mass ratio, the potential barrier at the cusp 
(or in an equivalent manner, the initial mass accretion rate), 
or the efficiency of the angular momentum transfer from the 
disc to the black hole. The extended number of models con- 
sidered here permits to draw some conclusive answers on the 
likelihood of the runaway instability of non-self-gravitating 
discs, as well as to pin-down the critical value of the ex- 
ponent of the angular momentum distribution separating 
stable and unstable models. 

As we showed in Paper II, and we further clarify in 
the present investigation, the runaway instability is most 
likely avoided in non rigidly rotating discs as those formed 
in the coalescence of a binary neutron star system or in 
the gravitational collapse of a massive star. This result is 
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in agreement with recent pertu rbative analysis of vert ically 
integrated discs performed by iRezzolla et alJ ll2003f) who 
found that the innermost regions of discs tend to behave as 
regions of evanescent-wave propagation, reducing the mass 
flux and, as a result, suppressing the runaway instability. If 
the growth of non-axisymmetric modes in the disc by the 
Papaloizou-Pringle instability is su ppressed by the accre- 
tion process itself, as suggested by Blae^ il987l) . systems 
consisting of a Kerr black hole surrounded by a high den- 
sity torus may then be long lived. The lifetime is probably 
controlled by the viscous timescale (a few seconds) rather 
than by the dynamical one. This may provide enough time 
for any plausible magneto-hydrodynamical process to effi- 
ciently transfer part of the energy reservoir of the system to 
a relativistic outflow. Therefore, the most favoured current 
models of GRBs (see IWooslevI i200 J) ) can indeed be based 
on such central engines. This is especially important in the 
context of the s o-called internal shock m odel for the prompt 
7-ray emission jRees fc Meszaro3ll994f) , where the observed 
lightcurve reflects the activity of the central engine, with no 
modification of the timescales other than the time dilation 
due to the redshift. In particular, the central engine has to 
survive for a duration at least comparable with the observed 
duration of the GRB. 

The paper is organized as follows: Section|5|describes in 
detail the relevant aspects of the theory of relativistic, sta- 
tionary, non-constant angular momentum thick discs neces- 
sary to build the initial state for our simulations. Part of 
this information is presented in three appendices at the end 
of the paper. Section |3 deals with a brief description of the 
numerical framework we use. Section |1| shows specific tests 
for the Kerr spacetime that the code has successfully passed. 
The description of the simulations and the discussion of the 
results are presented in Section |^ Finally Section |^ sum- 
marizes our conclusions and presents possible directions for 
future work. As we did in Papers I and II we use geometrized 
units (G = c = 1) unless explicitely stated. Usual cgs units 
are obtained by using the gravitational radius of the black 
hole Tg = GMbh/c^ = 1.5 x lO'^ (Mbh/Mq) cm as unit of 
length. Greek (Latin) indices run from to 3 (1 to 3), and 
we use a spacelike signature for the metric (- + + +). 



2.1 Gravitational field 

As mentioned before the self gravity of the disc is neglected. 
In Boyer-Lindquist {t, r, 6, 0) coordinates, the Kerr line ele- 
ment, ds^ = g^i.dx^dx'^ , reads 



2 A - sin^ 6,2 „ 2A/BHr sin^ 9 



ds = 



-dr - 2a- 



dtd(j> 



+ ^dr^ + Q^de^ + — Sin 



with the usual definitions: 



A 

2 

Q 



2,2 2 , 

r + a cos f 



(r + a ) — a A sin 9, 



(1) 



(2) 
(3) 
(4) 



where Mbh is the mass of the black hole and a is the black 
hole angular momentum per unit mass (Jbh/Mbh). We as- 
sume here that < a < A/bh so that prograde (respectively, 
retrograde) orbits correspond to positive (respectively, neg- 
ative) values of the angular momentum of the disc I (see 
below). The above metric, Eq. describes the spacetime 
exterior to a rotating and non-charged black hole. The met- 
ric has a coordinate singularity at the roots of the equation 
A = 0, which correspond to the horizons of a rotating black 
hole, r = r± = A^bh ± (Afln - a^Y''^- In the following, 
we use the notation rh = r+ and we call it "horizon of the 
black hole". Kerr black holes are also characterized by the 
presence of the ergosphere, a region inside which no static 
observers exist, whose boundary is given by 



r^{9) = A/bh + V ^-^BH - «^ cos- 



(5) 



Furthermore, the "distance to the rotation axis" is defined 
by 



Asm 9. 



(6) 



The 3+1 decomposition (see, e.g.. iMisner et all ((iS 

of this form of the metric leads to a spatial 3-metric 7ij with 
non-zero elements given by jrr ~ g'^/Ay^yee = g'^j'tr/itf, ~ 
Y^l (? SVC? 9. In addition, the azimuthal shift vector (i^ = gt^, 
is given by 



2aAfBHr'sin 9 



(7) 



2 EQUILIBRIUM INITIAL TORI 

The theory of stationary, relativistic thick discs (or tori) 
wit h a baryotropic equat ion of stat e (EoS) was derived 
by 'Kozlowsk i et al.l (|l97Sh (see also iFishbone fc Moncried 
(1976) ). Wc included in Paper I all relevant aspects of this 
theory necessary to construct the initial state for our simu- 
lations. In particular we described in detail the procedure to 
build equilibrium configurations of the system in the case of 
a Schwarzschild (non-rotating) black hole and a constant an- 
gular momentum disc. Here we extend this procedure to the 
more general case of a Kerr black hole and a non-constant 
angular momentum disc which follows a power-law distribu- 
tion with the radial distance. The notations and conventions 
we use are the same than in Paper I. 



and the lapse function is given by 



(8) 



For convenience we assume in the following that the 
initial mass of the black hole equals unity (therefore < 
a < 1). 



2.2 Angular momentum 

The angular momentum per unit inertial mass (hereafter 
"angular momentum") of a fluid element is related to the 
four-velocity by 



Utf, 
Ut 



(9) 



The corresponding angular velocity Q, is 
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n : 



gt4, + gtti 

g<t>4> + 9t4>l 



(10) 



The energy per unit inertial mass of a fluid element —ut can 
be obtained from the condition = Unu'^ = ~1: 



gttP + 'igt4,l + 04,4, 

This imposes a condition on the angular momentum: 

gtd^ + 2gt4,l + g^^, > , 

which is equivalent to 

I < l+{r, e)orl> l~{r, 6) for r < r^{e) 
l^, {r,e) <l<l+(r,e) for r > r^(e) 



(11) 



(12) 



(13) 



with 

~gtt 

A fluid element will be bound if 
equivalent to 



(14) 

'lit < 1, which is 



l+{r,9)< I 01 I >l-{r, 
l-{r,e)<l<l+{r, 



for r < ra{6) 
for r > re{9) 



with 



. ± UJy^l + gtt 

-gtt 



(15) 



(16) 



Notice that inside the ergosphere (r < re (6')), we have 
< 'cr < 'cr < ^'^^ outside the ergosphere (r > rc{9)), 
we have l^^ < < < l^^. On the surface of the ergosphere 
(r = re{6)), — l~. = —00, = (r^ + a^)/2a and = 
(r^ — rc + a^)/a. Then the angular momentum corresponding 
to the marginaUy bound case fuUfils condition 1131 . At the 
horizon (r = rh), we have = itr = 2rh/a. 



2.3 Keplerian orbits 

In the particular case of a Keplerian orbit in the equatorial 
plane, the angular momentum is given by 



K,cq 



= ±- 



(r- 2)7? ±a 



(17) 



Throughout the paper the subscript 'eq' means that the 
considered quantity is evaluated at the equatorial plane 
9 = n/2. In Eq. 11711 stands for a prograde (+) or 

a retrograde (-) orbit. A Keplerian orbit will be stable for 
d\lK,cci\/dr > 0. The radius of the last (marginally) stable 
orbit then corresponds to the case where dlK,aq/dr = 0, i.e. 
the minimum (respectively, maximum) of Zx.eq for a pro- 
grade orbit (respectively, retrograde orbit). This radius is 
given by 



■ 3 + y/{3~ Zi){3 + Z1 + 2Z2) 



(18) 



where ^'i = 1 + (1 - 0")^/=* [(1 + a^^"" + (1 - a)'/^] and 

Z'z — ^ Za? + Z\. A Keplerian orbit will be bound if /K,oq 
verifies equation 1)15^ . The radius of the last (marginally) 
bound orbit corresponds to the case where l^^^ = ^beq- 
This radius is given by 



2 =F a + 2\/l T a 



(19) 




0.4 0.6 
Black hole spin a 



Figure 1. Characteristic radii in the equatorial plane: the 

radius of the horizon r^, of the last marginally bound Keplerian 
orbit r^j^ and of the last marginally stable Keplerian orbit 
are plotted as a function of the spin a of the black hole. The 
radius r~ where the Keplerian angular momentum of a parti- 
cle with a retrograde orbit becomes infinite is also plotted. The 



schwarzschild limit (a = 0) corresponds to 



2, r 



and r. 



± 

ms 



mb 



6. The maximum black hole rotation limit (a = 1) 



corresponds to 
'"ins = 9 and ; 



1, 



2V2 ~ 5.83, 



(3 -I- \/5)/2 ~ 2.62. 



Notice that outside the horizon, the Keplerian angu- 
lar momentum becomes infinite for retrograde orbits at a 
radius r^ > 2 solution of (r — 2)-^ — a — 0. The charac- 



teristic radii rh. 



and r^r are plotted as a function 



of a in Figure All radii are given in units of the black 
hole gravitational radius. Correspondingly, the characteris- 
tic values of the angular momentum in the equatorial plane 
oq' ^ oq ^^'-^ 'cr.cq ^re plotted as a function of r in Figure|5| 
for a = (non rotating black hole), a = \/5/3 and a = 1 
(maximally rotating black hole). The intermediate value has 
been chosen to have rh = 1 + 2/3, where deviations from the 
Schwarzschild case start being more prominent. Notice that 



at the horizon 



b,cq 



2rh/a. 



2.4 Rotation law in the disc 



The distribution of angular momentum in the disc is as- 
sumed to follow a power-law in the equatorial plane: 



/cq(r) = Kr° 



(20) 



with Q > 0. Prograde (respectively, retrograde) discs corre- 
spond to K > (respectively, K < 0). From the rotation 
law in the equatorial plane, we can compute the value of 
the angular momentum at any location (r, 9) in the disc, 
by solving the equations of the constant I surfaces, i.e. the 
so-called von Zeipel's cylinders: 
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Figure 2. Characteristic values of the angular momentum in the equatorial plane: the critical value Zci,oq for which the energy 
per unit inertial mass of a fluid element —ut becomes infinite is plotted as a function of r in thick solid line and defines a forbidden 
region for the angular momentum which is indicated (grey area). The angular momentum for which the orbit is marginally bound 

is plotted in thin solid line. The region where orbits are not bound is also indicated (between thick and thin solid lines). The angular 
momentum of a Keplerian orbit is plotted in dashed line. The radius of the horizon rjj, the ergosphere r'c, the marginally bound 

Keplerian orbit and the marginally stable Keplerian orbit are indicated, as well as the radius r^i where the Keplerian angular 
momentum of a retrograde orbit becomes infinite. The interior of the black hole horizon has been shaded with crossing lines. 



l{r,e) = «oq(»"o), 



(21) 



where ro is the radius at which the von Zeipel's cylinder 
passing at {r, 6) intersects the equatorial plane. The struc- 
ture of these cylinders is presented in detail in Appendix fHl 
For practical purposes the value of l{r,6) = Zcq('"o) is ob- 
tained by solving numerically equation lB3t to find ro. 



2.5 Iso— pressure and iso— density surfaces 

The integral form of the relativistic Euler equation that gov- 
erns the dynamics of the gas flow can be written as 



W -Wi^ = \n{-ut) - ln(- 



-Ut 



ni 



(22) 



The subscript 'in' here refers to the inner edge of the disc in 
the equatorial plane. The "potential" W is deflned by 



(23) 



where p, p and h are the pressure, the density and the specific 
enthalpy of the gas (comoving frame). We limit ourselves 
here to the case of an isentropic fluid, where 



dp ^ h 

ph hin 



(24) 



In this case, the equipotentials W — const coincide with the 
equi-pressure p = const and the equi-density p = const sur- 
faces. Once the structure of the equipotentials W{r, 9) has 
been found, the thermodynamical quantities can be easily 
deduced through the condition h = /lin exp (Win — W) and 
the equation of state (EoS). In this paper, we will consider 
the case of polytropic fluids where 



p = up 
and 



(25) 



7-1 

so that 



P = 



p = 



- 1 (ftin e^'"-^ - 1) \ 



^-1 fone^.--^-l) \ 



(26) 

(27) 
(28) 



These expressions of h, p and p are valid only in the interior 
of the disc, which is defined hj W > Win. The total mass in 
the disc is given by the integral 



My: 



2-K 



g4><t> — gtti 



4- 2gt,pl + gttl' 



■ {ph + 2P) 



X 



p>0 



■ 6* ) sin 9 drd9. 



(29) 



In the models we consider in this paper, we have usually 
£ < 1 so that ph + 2P ~ p. 

According to these expressions, the construction of the 
equilibrium configuration of a thick disc reduces to the cal- 
culation of the function W{r,9). In the case of a constant 
angular momentum in the disc (see Paper I), the last term 
in Eq. 1221 vanishes and W{r,9) reduces to the first term, 
the estimation of which is trivial through Eq. Ijll^ once the 
distribution l(r,9) has been computed as described in the 
previous subsection. Here we consider more general distri- 
butions of I for which all terms have to be taken into account. 



2.6 Geometry of the equipotentials 

It is interesting to describe the structure of the equipoten- 
tials in the general case where a 7^ (the constant angu- 
lar momentum case a = has been fully described in Pa- 
per I). In the equatorial plane, the expression of Wcq(r) = 
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a = a = \/5/3 a = I 




Figure 3. Classification of the possible geometries of the equipotentials for three representative values of the black hole 
spin a. In the parameter space K-a of the rotation law lcq{r) = Kr°' in the equatorial plane, the three lines K = Kcv (thick solid line), 
K = KmB (thin solid line) and K = (dotted line) define three regions. In region 1, there is neither a cusp nor a centre; in region 

2, there are a cusp and a centre and the equipotential of the cusp is closed, which defines a disc; finally, in region 3, there are a cusp 
and a centre (the centre being rejected at +oo for a > 1/2) but the equipotential of the cusp is not closed, so that the disc is infinite. 
The dashed region in the plots is not considered in this paper as the derivative of the potential in the equatorial plane is not defined 
everywhere outside the horizon. 



W{r,n/2) can be written as 

f -l-oo 



H/cq(r) = - 



dr 



Wcq{r)dr 



1 



(30) 



where we have used the condition W{r,9) — > Oforr +oo 
to eliminate the radius rin and the corresponding value Win- 
From the definition of the von Zeipel's cylinders, the poten- 
tial at any given position is related to the potential in the 
equatorial plane by 



W(r,e) = W,q{ro)+ln 



-ut{r,e) 
-Ut{ro,Tv/2) 



(31) 



where ro is the radius where the von Zeipel's cylinder pass- 
ing at (r, 9) intersects the equatorial plane (ro is solution of 
Eq. (I2H ). The closed (respectively, open) equipotentials cor- 
respond to 14^ < (respectively, W > 0) and the marginal 
equipotential W = is closed at infinity. In order to un- 
derstand the geometry of the equipotentials, it is useful to 
study the derivative Wcq{r) = dWcq/dr which is given from 
equations ll^ . (ITTl and lHUIl by 



Weq(r) 



dr 



as 



j2 I r, 8gt4, j 



dr 



gttlcq + 'igttfilcq + §4,4, 



(32) 



Hereafter, the notation gtt, gt4>, §4,4, and zu means that these 
quantities are computed in the equatorial plane 9 — 7r/2. 
The function given by Eq. 1321 is studied in Appendix \C\ 
The condition for it to vanish is 



lcq{r) = ^K,cq('') , 



(33) 



which is not surprising: by the definition of a Keplerian orbit, 
when the angular momentum equals the Keplerian value, the 
gravitational forces are exactly compensated by the centrifu- 
gal forces and the pressure gradient vanishes; the pressure 
(and the density) are either minimum (surface of the disc) 



or maximum (centre of the disc) . Therefore the geometry of 
the equipotentials is mainly fixed by the distribution lcq{r). 
We now describe the geometrical structure of the equipoten- 
tials W = const for the particular choice of the distribution 
of angular momentum made in this paper (see Eq. 1201 ). 

2. 6. 1 Allowed values of the constant K 

We restrict this study to the case where the derivative of 
the potential Wcqij-) in the equatorial plane has only finite 
values outside the horizon. This corresponds to the fact that 
equation lcq{r) = iti.cqi'f') has no solution (see Appendix 1^. 
This is equivalent to 

\K\ < , 

where the marginal case K = Kcr is given by 



(34) 



^cr,cq(^cr) 



The radius is the solution of 

d In I /cr.cq I / \ 

— — (rcr = a . 

dlnr 



(35) 



(36) 



The slope of |^„,cql is an increasing function of r. For pro- 
grade orbits, the slope at the horizon is 

dln\l+^cq{r)\ _ { -00 0<a<l 



lim 

r — »r 



a = 1 



(37) 



and for retrograde orbits, the slope at the ergosphere is 

d\n\l^^^^^{r)\ _ 



lim 

7 — .2 



din r 



In both cases, the maximal slope equals 
d\n\lt,^q{r)\ _ 



lim 

r — » + oo 



dlnr 



(38) 



(39) 



Then, there are four cases: (i) for < a < 1/2 (sub- 
Keplerian case) and < a < 1, there is one unique root rcr 
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outside the horizon (respectively, outside the ergosphere) for 
prograde (respectively, retrograde) orbits except for a = 1 
and prograde orbits. In this last case, the root of equa- 
tion I36I I is inside the horizon and has to be replaced by 
r-cr = rh = 1 so that = 2; (ii) for 1/2 < a < 1 (Ke- 
plerian and super- Keplerian case) and < a < 1, there is 
one unique root rcr outside the horizon (respectively, out- 
side the ergosphere) for prograde (respectively, retrograde) 
orbits; (iii) for q = 1, this unique root is rejected at infinity 
so that rcr = +oo and Ka = 1 (respectively, —1) for pro- 
grade (respectively, retrograde) orbits; (iv) for a > 1, equa- 
tion ^eq('') = ^(!r,eq('') has always a solution. For this reason, 
we will not consider the a > 1 case in this paper. Notice 
that the case presented in Paper I corresponds to a = and 
a = 0. In this case r^r = 3 and Ka = 3\/3 ~ 5.20. 



2.6.2 Presence of a cusp and centre 

The geometry of the equipotentials is mainly determined by 
the roots of the equation lcti{r) = /j^ cq('')- condition to 

corresponding respectively 



have two roots rcusp and Tea 
to the cusp and the centre of the disc, is 

where the marginal case K — K^ia is given by 
(r-ms) 



The radius rms is the solution of 



dlnr 



(rms) = a. 



(40) 



(41) 



(42) 



The slope of |/jjj,q| is an increasing function of r. It equals 
at the radius of the last marginally stable orbit rms- The 
maximal slope equals 



lim JT"^ 

I — ►+CX) dlnr 



[r)\ 



(43) 



so that for < a < 1/2 (sub-Keplerian case) the root r^ns is 
unique and finite, for a = 1/2 (Keplerian case) the root rms 
goes towards infinity {Kms ~^ ±1), and for a > 1/2 (super- 
Keplerian case) there is no root: Eq. I|33^ has always only 
one root, the radius of the cusp rcusp. Notice that for a — 
(constant angular momentum, see Paper I), the root rms 
coincides with the the radius of the last marginally stable 
orbit and Kms = fe,eq(T'ms) in this case. 



2.6.3 Closed equipotential at the cusp 

In the case where there is a cusp at rcusp the corresponding 
equipotential will be closed only if 



W^cq(?-cusp) < 0. 

This is equivalent to 

'^^t.cq (rcusp ) ^ e ^ ^ . 



(44) 



(45) 



In case of equality, the equipotential of the cusp is closed at 
infinity (VKcusp ~ 0). The quantity Fcq{r) is defined as 



1 ^cq^cq dv 



(46) 



gttl'iq 



gt4, + gttl 



cq dlcQ 

^ -dr. 

- 2(/i^<cq + dr 



(47) 



As already mentioned, in the case where a = the integral 
Jcq(r) vanishes and this condition is equivalent to the condi- 
tion for a fluid element in orbit at rcusp with angular momen- 
tum ;cq(rcusp) to be bound (Eq. (O), i.e. \K\ < |^K,cq(''mb)l- 
In the more general case it can be shown that the condition 
is equivalent to 



\K\ < \Kn,b\, 

where Kmb is the solution of 



W^cq(rc 



for = ATmb. 



(48) 



(49) 



For < Q < 1/2, we have \Kcr\ > |A'mb| > |ifms| > 1 (with 
-K^mb = Kms = ±1 for a = 1/2): the cusp and the centre are 
respectively a local maximum and minimum of the potential 
Wcq{r). For a > 1/2, the centre is rejected at infinity where 
M4q(r) and the cusp is still the maximum of VKcq(r) 
and then must have a positive potential. This means that 
^mb = in this case. 

Figure 13 shows the critical coefficients Kct, Kmh and 
Kms as a function of a for pro- and retrograde discs and for 
three different values of the spin of the black hole: a = 
(non-rotating), a = y/E/3, and a = 1 (extreme rotation). 
This allows to classify the possible geometries of the equipo- 
tentials, which is presented in TableQ This table is the gen- 
eralized version of Table 2 in Paper I. The lines appearing in 
Fig. 13 separate different equipotential geometries. In region 
1, there is neither a cusp nor a centre. Correspondingly, in 
region 2, there are cusp and centre and the equipotential 
of the cusp is closed, which defines a disc. Finally, in re- 
gion 3, there are cusp and centre but the equipotential of 
the cusp is not closed, so that the disc is infinite. The cor- 
responding cases are illustrated in Figures 2] (for a = 0),|K| 
(for a — \/5/3), and El (for a — 1). We notice that in order 
to limit the number of figures to a reasonable amount we 
do not include the case of retrograde discs/rotating black 
holes. Figs. I4I6I show the various equipotential geometries 
that may arise depending on the value of K, and this is done 
for the sub-Keplerian (a < 1/2), Keplerian {a = 1/2) and 
super- Keplerian {a > 1/2) cases. In all these figures, the left 
panels show the radial distribution of the angular momen- 
tum in the equatorial plane and the corresponding potential. 
In addition, the right panels show the equipotentials in the 
X — y plane. In the case a > and a retrograde disc (i.e. 
if < 0), the evolution of the geometry of the equipotentials 
with 17^1 is essentially identical and we do not include it 
here. The main difference is that the cusp - when it exists - 
is at a larger distance from the horizon, as well as the centre 
of the disc in case (2). The interesting case for the study of 
the runaway instability corresponds to the situation when 
there are both a cusp and a centre, and where the equipo- 
tential of the cusp is closed, which means < a < 1/2 and 
\Kms\ < \K\ < l-K^mbl. Notice that the size of the disc evolves 
from zero {K — K^s) to infinite {K — K^b), which allows 
to fix its mass by adjusting the constant K. This case of 
interest is illustrated in the third row of the first page of 
Figs. 14161 In this case, we define the potential barrier at the 
surface of the disc by 



AWin = Win - Wcusp. 



(50) 
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Table 1. All possible equipotential configurations for a Kerr black hole and a disc with an angular momentum distribution in the 
equatorial plane following a power- law le<i{r) = Kr". 

Prograde rotation : K > 



sub-Keplerian case < a < 1/2 



Angular momentum 


fcusp 


Wcusp 


''centre 


^centre 


Comments 


0<K < Kms 
K = K^s 
Kms <K < 
K = K^b 

Kmb <K <Kcr 
K = Kcr 


''cusp — '^ms 
^mb ^cusp '^ms 
T'cusp = fjnh 
T'cusp *C fjah 


< 

< 


> 


''centre — ''cusp 
''centre ^ fras 
''centre ^ ''ms 
''centre ^ ''ms 
''centre ^ ''ms 


o o o o o 

V V V V V 


No cusp. No centre. No closed disc. 
Cusp = centre. Disc reduced to one point. 
Cusp. Centre. Disc closed. 
Cusp. Centre. Disc closed at infinity. 
Cusp. Centre. Disc infinite*. 
No cusp. Centre. Disc infinite*. 








Keplerian case a = 


1/2 




Angular momentum 


^cusp 


Wcusp 


''centre 


^centre 


Comments 


<K <1 
K = 1 

KK <Kcr 
K = Kcr 


+00 
rcusp <+oo 



> 


H-CX) 
+CX) 
+00 







No cusp. No centre. No closed disc. 
Cusp and centre at infinity. No closed disc. 

Cusp. Center at infinity. No closed disc. 
No cusp. Center at infinity. No closed disc. 


super-Keplerian case 1/2 < a < 1 


Angular momentum 


Tcusp 


Wcusp 


^centre 


VKcentre 


Comments 


0<K <Kcr 

K = KcT 


rcusp < +00 


> 


+ 00 
+CX) 






Cusp. Center at infinity. No closed disc. 
No cusp. Center at infinity. No closed disc. 






Retrograde rotation 


: K <Q 








sub-Keplerian case < 


a < 1/2 




Angular momentum 


Vcusp 


l^cusp 


''centre 


Wcentre 


Comments 


Jfms < -fST < 
K = K^s 
Kmh <K < Kms 
K = K^b 

Kcr<K< 


''cusp — t^ms 
fmh ''cusp <^ ^ms 
Tcusp = fmh 
''cusp *^ ''mb 


< 

< 


> 


''centre — ''cusp 
''centre ^ ''ms 
''centre ^ ''ms 
''centre ^ ''ms 


o o o o 

V V V V 


No cusp. No centre. No closed disc. 
Cusp = centre. Disc reduced to one point. 
Cusp. Centre. Disc closed. 
Cusp. Centre. Disc closed at infinity. 
Cusp. Centre. Disc infinite*. 








Keplerian case o; = 


1/2 




Angular momentum 


Tcusp 


Wcusp 


''centre 


VKcentre 


Comments 


-1< A- < 
K = -1 

Kcr<K< -1 


+00 
rcusp < +00 



> 


H-CX) 
+CX) 






No cusp. No centre. No closed disc. 
Cusp and centre at infinity. No closed disc. 
Cusp. Center at infinity. No closed disc. 


super-Keplerian case 1/2 < a < 1 


Angular momentum 


''cusp 


VKcusp 


''centre 


Weentre 


Comments 


Kcr<K<0 


rcusp < +0O 


> 


+ 00 





Cusp. Center at infinity. No closed disc. 



* Some closed equipotentials are still present around the centre. 



There are three possible cases to consider: (i) for ATVin < 0, cusp; (iii) the marginal case, AWin = 0, which corresponds 
the disc is inside its Roche lobe (defined as the equipotential to the case where the disc is exactly filling its Roche lobe, 
of the cusp) and there is no possible mass transfer from the 
disc to the black hole; (ii) for AWin > 0, the disc is overfilling 
its Roche lobe and mass transfer is possible through the 
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Schwarzschild black hole (a = 0) 



Sub-Keplerian case (here with a = 1/4) 
Case (1): K < K^^s (here K = ±1.5 ; no cusp ; no centre) 




Figure 4. Geometry of the equipotentials: Schwarzschild black hole (a = 0) and pro- or retrograde disc. Sub-Keplerian 
case. For each case, the angular momentum Icqii") and the potential Wcq{r) in the equatorial plane are plotted in the left panel, and the 
equipotentials in the x = r sin 6-y = r cos 6 plane in the right panel. Left: the critical (respectively, Keplerian) angular momentum i^.cq (*") 
(respectively cq^*^')) plotted in thin solid (respectively, dotted) line. The intersections of Zeq(''') and lK.cq{i") (if any) correspond to 
the cusp or to the centre, which are respectively a local maximum or minimum of Weq(r). Right: the interior of the black hole horizon is 
filled in black. The critical von Zeipel cylinder is indicated with a dashed line. The equipotentials outside the critical cylinder are plotted 
in dotted line and are extrapolated using a constant angular momentum in this region, equal to the value on the critical cylinder. When 
present, the equipotential of the cusp is plotted in thick line and the centre is indicated with a big dot. 
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— Schwarzschild black hole (a = 0) — 

Sub-Keplerian case (here with a = 1/4) 
Critical case: K = K^i, (here K^\, ~ ±2.35; rcusp ^ 5.96 ; rcentre — 20.1) 




10 102 
Radius r / 

Critical 



10 15 

X / 

±3.90; no cusp ; rcentre — 223) 




10 10^ 
Radius r / 



Figure 4 - continued Geometry of the equipotentials: Schwarzschild black hole (a = 0) and pro- or retrograde disc. 
Sub-Keplerian case. In the critical case K = Xmax (third line), the thick line corresponds to the critical equipotential W = +oo. 
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— Schwarzschild black hole (a = 0) — 

Keplerian case (a = 1/2) 
Case (1): K < K^^s = ^t^mb = 1 (here K = ±0.5; no cusp ; no centre) 



5 10 




10 102 
Radius r / 

Case (3): 1 = Kms = -?^mb < K < A'max (here K = ±1.5 ; rcusp = 6; centre at oo) 



5 10 




10 102 
Radius r / 




Figure 4 - continued Geometry of the equipotentials: Schwarzschild black hole (a = 0) and pro- or retrograde disc. 
Keplerian and super-Keplerian cases. 
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- Mildly rotating Kerr black hole [a = \fE/Z) - 

Sub-Keplerian case (here with a = 1/4) 
Case (1): K < -ftTms (here K = 1.5; no cusp; no centre) 




1 10 10^ 

Radius r / 

Critical case: K = A'ms (here Kms — 1.99; rcusp = r-contrc ~ 5.45) 




10 10= 
Radius r / 

Case (2): Kms <K < K^^ (here K = 2.02; rcusp ^ 4.05; rccntrc ^ 7.74) 




Figure 5. Geometry of the equipotentials: mildly rotating Kerr black hole (a = y/Z/Z) and prograde disc. Sub-Keplerian 
case. All line styles and notations are identical to those used in Fig. |1] 



The runaway instability of thick discs around black holes. II. 13 

- Mildly rotating Kerr black hole (a = \/5/3) — 

Sub-Keplerian case (here with a = 1/4) 
Critical case: K = K^i, (here K^\, ~ 2.09; rcusp — 3.29 ; r-centre — 10.8) 




1 10 10= 10= 5 10 15 20 25 

Radius r / x / r^ 



Case (3): K^h < K < JCmax (here K = 2. ~ 2.76 ; rcentre — 15.6) 




Radius r / r^ x / r^ 

Figure 5 - continued Geometry of the equipotentials: mildly rotating Kerr black hole (a = \/5/3) and prograde disc. 
Sub-Keplerian case. 
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- Mildly rotating Kerr black hole (a = \/5/3) — 

Keplerian case (a = 1/2) 
Case (1): K < /sTms = ^S^mb = 1 (here K = 0.5; no cusp; no centre) 




-0.2 



Case (3): 1 = Kms = i^mb < K < Kmax (here K = 3.54; rcusp = 6; centre at oo) 




- 



10 10^ 
Radius r / 



20 25 




Figure 5 - continued Geometry of the equipotentials: mildly rotating Kerr black hole (a = VE/S) and prograde disc. 
Keplerian and super-Keplerian cases. 
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— Extreme Kerr black hole (a = 1) — 

Sub-Keplerian case (here with a = 1/4) 
Case (1): K < /■Cms (here K = 1.5; no cusp; no centre) 




Figure 6. Geometry of the equipotentials: extreme Kerr black hole (a = 1) and prograde disc. Sub-Keplerian case. All 

line styles and notations are identical to those used in Fig. |1] Note the scale change in the plots on the right column as compared to 
Figs. El and El 
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c 10 



I 1 





-0.02 



--0.04 
o 

-0.06 



— Extreme Kerr black hole [a = 1) — 

Sub-Keplerian case (here with a = 1/4) 
Critical case: K = K^i, (here — 1-86; rcusp — 1-52 ; r-centre — 4.87) 

6 




I I I I I III 1 I I I I III 1 M il - 




1 10 10= 

Radius r / r 



Case (3): K^^ <K < K^^ (here K = 1.87; rcusp ^ 1.47 ) '^centre — 5.07) 



c 10 





10 10= 8 4 6 8 10 

Radius r / x / 

Critical case: K = Jfmax (here Kmax = 2.; no cusp; rcentre — 9.18) 



c 10 ^ 




O'-0.02 



10 10= 
Radius r / r^ 




Figure 6 — continued Geometry of the equipotentials: extreme Kerr black hole (a = 1) and prograde disc. Sub-Keplerian 
case. 



The runaway instability of thick discs around black holes. II. 



— Extreme Kerr black hole (a = 1) — 

Keplerian case (a = 1/2) 
Case (1): K < JsTms = ^mb = 1 (here K = 0.5; no cusp; no centre) 

6 




Case (3): 1 = K^s = K^^, < K < -Kmax (here K = 3.54; rcusp — 2.17; centre at oo) 

6 



3 

C 10 




0.4 
^ 0.2 

I -0.2 
-0.4 



ft \ — I I I III 1 — h 




-4 



10 IQS 
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Super-Keplerian case (here with a = 3/4) 
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Case (3): K < ii'max (here K = 1; rcusp — 3.32; centre at oo) 




10 IQS 
Radius r / 




Figure 6 - continued Geometry of the equipotentials: extreme Kerr black hole (a = 1) and prograde disc. Keplerian 
super-Keplerian cases. 
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2.7 Method of construction 

A configuration of a thiclc disc orbiting a black liole is defined 
by eiglit parameters: (i) two parameters describing the black 
hole, the mass Mbh and the spin < o/Mbh < 1; (ii) four 
parameters for the disc: the mass Md, the slope < a < 1/2 
and the sign (prograde or retrograde rotation) of the angular 
momentum in the equatorial plane, and the potential barrier 
at the surface A Win; (iii) two parameters for the EoS: the 
adiabatic index 7 and the polytropic constant k. Once this 
set of parameters has been specified, we adopt the following 
procedure to build the disc: 

(i) Gravitational field: from a, we compute the Kerr met- 
ric coefficients appearing in Eq. Q. 

(ii) Limits for the angular momentum: from a and the 
sign of K, we compute the two limits ifms and K-^h defined 
by equations II4H and 1491 1. Notice that /sTmb is computed 
by solving Eq. 1491 with a bisection procedure, where each 
iteration implies the computation of rcusp and Wcq(rcusp). 

(iii) Initial value of K: |/('ms| < \K\ < li^mbl- 

(iv) Potential in the equatorial plane: the function Wcqir) 
is computed from Eq. 13011 . 

(v) Characteristic radii: we compute the radius of the 
cusp and the centre rcusp and rcontro by solving Eq. 13311 . 
and the inner radius rin from AVKin. 

(vi) Energy per unit inertial mass: the function ~-ut{r,0) 
is computed everywhere from Eq. IIH . where the angular 
momentum l{r, 9) is computed from the procedure described 
in Subsection 12.41 and Appendix lEl 

(vii) Potential: the function W{r,6) is computed from 
Weq{r) and —Ut(r,B) using Eq. ll^ . 

(viii) Enthalpy, density and pressure: h{r,6), p{r,6) and 
p{r,6) are computed from equations 126L 1271 and 12811 in 
the region where W < Win (interior of the disc) and are 
forced to h = 1, p = 0, p — elsewhere (vacuum). 

(ix) Mass: Md is computed from the integral 1291 . 

(x) From the difference between Mo and its expected 
value, the coefficient K is adjusted (secant method) and the 
iteration starts again at step (iv). This uses the fact that 
Md/A/bh as a function of \K\ is strictly increasing from 
for K = Kms to +00 for K = K^ih- This procedure is 
stopped once the accuracy obtained for the value of Mu is 
good enough. 



3 NUMERICAL FRAMEWORK 

The axisymmetric numerical code we use to evolve in time 
the equilibrium configurations described in the previous sec- 
tion is the same we used in Paper I. The code was orig- 
inally written to account for the Kerr metric terms. The 
Schwarzschild black hole case analysed in Paper I was hence 
studied by setting to zero the black hole angular momen- 
tum parameter. In our code the general relativistic hydro- 
dynamics equations (see Paper I for the especific expre- 
sions in Kerr spacetime) are numerically solved using an 
upwind high-resolution shock-capturing scheme based upon 
the HLLE Riemann solver. The code is second order accu- 
rate in both space and time due to the use of a piecewise 
linear, cell-reconstruction algorithm and a two-step, conser- 
vative Runge-Kutta time update. 

The form of the Kerr metric we employ is the one 



given by standard Boyer-Lindquist (spherical) coordinates 
{t,r,9,(f)) (cf. Eq. (Q). Axisymmetry with respect to the 
black hole rotation axis implies that all (f) derivatives van- 
ish. The computational grid has 400 x 100 zones in the radial 
and polar direction, respectively, and it is logarithmically 
spaced in the radial direction to ensure the finest resolution 
at smaller radii. The typical width of the innermost cell, 
where we have the highest resolution, is Ar ~ 1.99 x 10~^. 
The innermost radius is located at Tmin = 2.12, sufficiently 
close to the horizon of the initial Schwarzschild black hole 
to avoid numerical inaccuracies introduced by the coordi- 
nate singularity of the coordinate system we use. It is clear 
that the use of " horizon-adapted" coordinate systems as the 
ones adopted bv lFont et alJin'999ll . regular at the black hole 
horizon, would suppress all numerical difficulties associated 
with the position of the innermost radius of the grid. This 
approach, however, has not been pursued in the present in- 
vestigation but will be considered in the future. We postpone 
the explanation of the procedure we employ to increase the 
innermost radius in simulations where the size of the black 
hole horizon is growing to Section [5.21 below. On the other 
hand, the location of the maximum radius rmax depends on 
the model. For the stationary test models presented in Sec- 
tion 2] rinax = 35. Correspondingly, for the simulations of 
the runaway instability, the radial grid extends to a suffi- 
ciently large distance in order to ensure that the whole disc 
is included within the computational domain. The partic- 
ular values of rmax that we choose are reported in Table |5| 
below. Finally, as we did in Paper I, we use a finer grid in the 
angular direction within the torus and a much coarser grid 
outside. The angular zones are distributed according to the 
same law presented in Paper I. The computational domain 
in the angular direction extends from to vr. 

Following the same approach than in Papers I and II we 
restrict here to adiabatic hydrodynamical simulations which 
reduces the number of equations to solve to four, the rela- 
tivistic counterparts of the continuity and Euler equations. 
This permits to obtain the specific internal energy through 
a polytropic EoS, p = i^p^ , i.e. e = t^p'''^ ■ The bound- 
ary conditions for the (primitive) hydrodynamical quantities 
also coincide with those discussed in Paper I. Furthermore, 
the procedure to recover those variables from the evolved 
quantities (the relativistic densities of mass and momenta) 
is the same of Paper I. We refer the reader to that work for 
further details. 

Finally, it is worth commenting on the simple yet effec- 
tive procedure we adopt in order to evolve with our hydro- 
dynamical code the "vacuum" zones which lie outside the 
disc. Before constructing the initial torus we build a back- 
ground, spherical 'dust' accretion solution of sufficiently low 
density so that its presence does not affect the dynamics of 
the disc. This solution corresponds to a radially accreting 
flow with negligible pressure (geodesic motion) . We even se- 
lect among all possible solutions the marginally bound case 
where ut = —1. Our 'dust' solution is then determined by 
the following equations, in which there is only one free pa- 
rameter: Ur = — \J — p = —j=- The pressure and 
the specific internal energy are computed from the density 
using a polytropic EoS. In our approach we adjust ifdust so 
that the maximum density in the background spherical so- 
lution is 5 X 10~® times the maximum density at the centre 



The runaway instability of thick discs around black holes. II. 



19 




10-2 0.1 1 10 10= 

t (ms) 

Figure 7. Time evolution of the mass flux for a stationary model 
with a Kerr black hole of spin q/Mbh =0.9 and a disc of con- 
stant angular momentum I = 2.6088. Models (a), (b), (c) and (d) 
correspond respectively to AVFin/c^ = 0.04, 0.08, 0.16 and 0.32. 
The time is given in ms and the mass flux in units of the Edding- 
ton mass flux, both assuming that the mass of the black hole is 
1 Mq . After about 1ms, the mass flux tends asymptotically to the 
stationary value. These test simulations are made with 7 = 4/3 
and re = 1.5 X 10-^" cgs. 

of the disc. By doing this we have checked that the rest- 
mass present in our background solution is always negligible 
compared to the mass of the disc and that the associated 
mass flux corresponding to this spherical accretion is also 
negligible compared to the mass flux from the disc. 



4 CODE TESTS 

As we did in Paper I for the case of a Schwarzschild black 
hole, we have first tested the code by considering time- 
dependent simulations of stationary constant angular mo- 
mentum discs around a rotating black hole. The aim of these 
simulations is to check whether the code is able to keep those 
tori in equilibrium during a sufficiently long period of time 
(much larger than the rotation period of the discs). In order 
to do so we ha ve considered the same four sta t ionary mod- 
els analysed bv lleumenshchev fc Beloborodovl il997l) . for a 
1 Mq rotating black hole with o/Mbh = 0.9. These four 
models are characterized by a constant angular momentum 
I = 2.6088, a value in between the marginally stable and 
marginally bound orbits, and an increasingly large value of 
the energy gap at the cusp, AVKin = 0.04, 0.08, 0.16, and 
0.32. The adiabatic index of the EoS is 7 = 4/3 and the 
corresponding polytropic constant k — 1.5 x 10^" cgs. Fur- 
thermore, the mass and the spin of the black hole are kept 
constant throughout these test evolutions. 

A conclusive proof of the ability of the code in keeping 
the stationarity of the various equilibrium models considered 
is provided by Figure |7| This figure shows the evolution of 
the mass accretion rate (normalized to the Eddington value), 
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Figure 8. Mass flux as a function of the potential barrier 
AWin/c^ ■ For each of the models of Fig.|3 the asymptotic station- 
ary value of the mass flux is plotted as a function of the potential 

barrier (big dots). The analytic dependence is rh oc AW-''^^ ■ 
The corresponding slope is 4 for 7 = 4/3 and is indicated with 
the straight line. 

wi/mEdd, as a function of time (in ms). The mass flux (flux 
of rest-mass density) is computed at the innermost radial 
point according to 

/•TT 

m = 27r / y/^Dv^'de, (51) 
Jo 

where D — pT, T being the Lorentz factor. The volume 
element for the Kerr metric appearing in the above equa- 
tion is given by y/—g = g^sind. The Eddington mass flux 
^Edd = ~ 1.4 X 10" g/s is computed for 

Mbh = 1 Mq. Rescaling of m for different polytropic con- 
stant K and black hole mass Mbh is given by 

(m/mEdd)i _ ^fea 1 

(m/mEdd)2 \ii2J Mbh 2 

where 7 = 4/3 is assumed. Fig. Q shows that after a tran- 
sient initial phase of about 1ms the mass flux rapidly tends, 
asymptotically, to a constant value. The offset observed dur- 
ing the initial phase corresponds to the spherical accretion 
mass flux associated with the particular background solution 
we use outside the torus. 

Correspondingly, Figure |H| shows the dependence of the 
mass fiux with the energy gap AVFin for the four stationary 
models we consider. The values selected for the mass ffux 
in each model are the asymptotic ones read off from Fig. |7| 
This plot allows us to check if the c ode is able to reproduc e 
the analytic dependence derived bv lKozlowski et al.l lll978l) 

m oc AW.J^ . (53) 

For a 7 = 4/3 polytrope the expected slope is 4 (dashed 
line). As the figure shows our results are in good agreement 
with this analytic prediction as well as with the numerical 
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results obtained bv llgumenshchev fc Beloborodovl (^O^) for 
the same models. 



5 RESULTS 

We turn now to describe our simulations. The main objec- 
tive of this section is to demonstrate the suppression of the 
runaway instability in thick discs around black holes due 
to the stabilizing effect of an increasing angular momentum 
distribution. 



5.1 The physical origin of the stabilizing effect 

As already mentioned in the introduction the runaway in- 
stability of thick d i scs ar ound black holes was discovered by 
lAbramowicz et all (119831) with a stationary model which em- 
ployed Newtonian gravity and a pseudo-potential to mimic 
the relativistic gravitational field of the black hole. In Paper 
I we summarized all subsequent studies which followed this 
first work, as well as the current understanding of the phys- 
ical origin of the instability. The results presented in Paper 
I were the first confirmation by hydrodynamical simulations 
in general relativity that the runaway instability is indeed 
occuring, at least for non self-graviting constant angular mo- 
mentum discs. Our study also provided the first confirma- 
tion that this instability was extremely violent, leading to 
the destruction of the disc on a dynamical timescale, shorter 
than 1 s for realistic par ameters. More recent ly these results 
have been confirmed bv lZanotti et alJ i2003l) with indepen- 
dent hydrodynamical relativistic simulations. These authors 
have also shown that the instability of constant angular mo- 
mentum discs is a robust feature, independent of the way 
accretion is induced. 

The main goal of the present study is now to check - 
using again hydrodynamical relativistic simulations - the 
stabilizing effect of a non-constant angular momentum dis- 
tribution in the disc. In the coUapsar model of GRBs, the 
simple constraint given by the Rayletgh's criterion imposes 
such a distribution both in the initial massive star before 
the collapse and in the final thick disc. The constant an- 
gular momentum case appears as a rather unlikely limiting 
case. Positive slope of the angular momentum in the disc 
are indeed found in numerical simulations of the gravita- 
tional collapse of a massive star, as well as in simulations 
of the coalescence of two compact objects. The stability 
of non- const ant angular momentum discs was first discov- 
ered bv lDaigne fc MochkovitchI 1^93) in pseudo-Newtonian 
gravity, and confi r med with a relativistic calculation by 
lAbramowicz et all (Il99^ . both studies being based on sta- 
tionary models. The physical origin of this stabilizing effect 
is simple: during the accretion process, the material at the 
inner edge of the disc is replaced by new material with a 
higher angular momentum. Therefore, the centrifugal bar- 
rier which acts ag ainst accretion becomes much more effi- 
cient. As shown bv lDaigne' fc MochkovitchI ((l993) the differ- 
ence between discs with constant and non-constant angular 
momentum is the location of the cusp (the Lagrange point 
Li) with respect to the inner edge of the disc after mass 
transfer. In the first case, the inner edge of the disc moves 
faster than the cusp towards the black hole horizon (favour- 



ing the instability) whereas the cusp moves faster towards 
the horizon in ^the latter (reg u lating the accretion). 



presented a perturba- 



Recently, Rezz olla et al 



five study of vertically-integrated discs around black holes. 
They found that when acoustic waves are considered, the in- 
ner region of the torus are of evanescent type for discs with 
non constant angular momentum, while such an evanescent 
region does not exist in the constant angular momentum 
case. This leads to a complementary explanation of the sup- 
pression of the runaway instability in non-constant angular 
momentum discs: the evanescent region in the inner disc 
prevents the inward propagation of perturbations and acts 
against the accretion of mass on to the black hole. 



5.2 The sequence of stationary metrics 
approximat ion 

The fiux of mass and angular momentum from the disc to 
the black hole increases the mass and the spin of the latter, 
therefore changing the metric potentials of the underlying 
gravitational field where the fluid evolves. Such increase is, 
in fact, the fundamental process which triggers the runaway 
instability of the disc. Thus, it needs to be included in some 
way in the hydrodynamical simulations. 

Strictly speaking, in order to properly take into account 
the dynamical evolution of the gravitational fleld, one should 
solve the coupled system of Einstein and hydrodynamics 
equations for a self-gravitating matter source. Numerical 
relativity codes capable of evolving black hole spacetimes 
with perfect fluid matter are onl y becoming available re- 
cently, both in axisymmetry (fBrandt et al. 2000; Shibat^ 
__05 ) as in full 3D (p hibata & UrvO 200(1 IShibata et all 
20od:lAlcubierre et aLjlToOQ : .Font et al...200lll . These codes, 
however, have not yet been applied to investigate astro- 
physical aspects of accretion on to black holes, in particular 
the runaway instability. When trying to simulate dynamical 
black hole spacetimes there are some issues which make this 
a challenging problem, among those the affordable resolution 
and the integration times required to study the development 
of the runaway instability, perhaps far too demanding for 
such relativistic codes which incorporate self-gravity. One of 
the most serious problem is, no doubt, the numerical han- 
dling of the physical singularity hidden inside the black hole 
horizon, which prevents long-term evolutions. Recently, suc- 
cessful, long-term simulations of (vacuum) black hole space- 
times have been achieved by excising part of the computa - 
tional grid within the event horizon (see e.g. lYo et al.l(l2002h 
and references therein). This is a promising line of research 
which could provide a way to investigate the runaway insta- 
bility in full general relativity in the near future. 

While these obstacles are paved away we have adopted 
a simplified and pragmatic approach to the problem. In our 
procedure the spacetime metric is approximated at each 
time step by a stationary exact black hole metric of varying 
mass and angular momentum. As we mentioned in Paper I 
the mass Mbh of the black hole, necessary to compute the 
metric coefficients, is increased at each time step At from 
time i" to time according to: 

M^+' = MSu + Atm" , (54) 
where the mass fiux at the inner radius of the grid is evalu- 
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ated by the equations 



dm 
U 



dm 

'de 



(55) 



(56) 



Similarly, the angular momentum Jbh of the black hole is 
increased at each time step according to: 



(57) 



where the angular momentum flux at the inner radius is 
computed by the equation 



3 = 1 



dm 



(58) 



The angular momentum is related to the (;A-component of 
the 3- velocity used in the hydrodynamics code by 



I = 



(59) 



In addition, the parameter rj appearing in Eq. 15811 controls 
the efficiency of the transfer of angular momentum from 
the disc to the black hole. Since in realistic discs angular 
momentum is transported outwards by dissipative processes 
or removed from the system by gravitational radiation, we 
adopt in most models a conservative value for the efficiency 
of the angular momentum transfer, assuming that only 20% 
of the angular momentum of the accreted material is trans- 
ferred to the black hole. 

An important numerical aspect is the evolution of the 
radial grid in simulations where the increase of the mass and 
the spin of the black hole is taken into account. The radius 
of the horizon Vh = Mbh + \/ -^'^bh ~ moves according to 



1 



Mbh 



2 

BH 



dMsH 



2 

BH 



-.da . (60) 



Following the procedure which has just been described, the 
variation of a is related to the variation of mass via the value 
lin of the specific angular momentum of the material passing 
through the inner limit of the grid and the efficiency ry of its 
transfer to the black hole: 

da _ dJsH _ dMBH _ / rjlm , \ dMaH 
a 



Jbh 



Mb 



Mb 



(61) 



By combining this expression with the previous equation 
and after some algebra we can express the variation of the 
radius of the horizon as 



drh 
rh 



^ 2 



Mbh ) l^T (^h) 

dMBH 



( Mbh ) 



Mbh 



(62) 



where l^^(r\^ — 2MbhT\iI a is the maximum value of the 
angular momentum at the horizon. This v ariation is positive 
for nhr^/ltirh) < l-l/2v'l-(a/A/BH)2, a condition which 
is always fulfilled in our simulations. Notice that for a = 
(Schwarzschild black hole), this condition reads r/hn < +oo 
and for a — Mbh (maximally rotating black hole), it reads 



rjiin < /^(rh). Therefore, in all of our simulations we observe 
the expected behaviour: the horizon moves outwards. Then, 
as we did in Paper I, we follow a simple procedure to let the 
radial grid evolve with time. We just increase the index imin 
of the innermost cell when necessary so that the condition 
^imin > is always respected. It is important to notice that 
for all disc-to-hole mass ratios considered in our simulations, 
the black hole mass and spin increase very slowly during the 
evolution. This implies that the metric coefficients at any 
time step differ very little from the final values which would 
correspond to an exact Kerr black hole of larger mass and 
spin but with no matter around. 

5.3 Initial state 

As shown in Section |2l a given initial state of the black hole 
plus non-constant angular momentum disc system is deter- 
mined by eight parameters: the mass and the spin of the 
black hole, the disc-to-hole mass ratio, the slope (and sign) 
of the angular momentum in the equatorial plane of the disc, 
the potential barrier at the disc surface, and the two parame- 
ters defining the polytropic EoS. In all simulations presented 
in this paper, the black hole is initialy non-rotating (a — 0) 
and therefore, the initial sign of the angular momentum of 
the disc is chosen to be positive (prograde disc) as the retro- 
grade case is strictly identical. The case of initially rotating 
black holes and pro- or retrograde discs will be presented in 
a forthcoming paper. As already mentioned in Paper I, the 
computing time needed for one hydrodynamical simulation 
is too large to allow for a complete exploration of the domain 
defined by the remaining six parameters. For this reason we 
focus on those models which are expected to be found in the 
central engine of GRBs. 

5.3.1 Black hole mass and disc-to-hole mass ratio 

In the two most discussed scenarios the central engine of 
GRBs is formed either after the coalescence of two compact 
objects or after the gravitational collapse of a massive star. 
As in Paper I, we then fix the mass of the black hole to 
be 2.5 M0, which is a reasonable value according to the 
results of various numerical simulations teuf fert fc_^ 
.199a: .Kluzniak fc Lee _^98; MacFadvcn & W ooslev 
Ishibata fc Urvul2000l: l 



Alov et al.l2000il:IShibata et al 



; to the 
IB; 



In addition we only consider two possible values for the 
disc-to-hole mass ratio, either 1 or 0.1. The second value is 
clearly more realistic according to the aforementioned sim- 
ulations (even some lower values ~ 10~^ — 10~^ are found 
in some situations). This small mass ratio makes our initial 
assumption of neglecting the self-gravity of the disc less of a 
concern. However, as the reservoir for the accretion process 
is larger in the case where Mo /Mbh = 1, the time neces- 
sary to completely destroy the disc in the case of a stable 
accreting regime is much longer than the timescale of the 
runaway instability observed in unstable models. Therefore, 
our choice of this high mass ratio is justified as it allows 
a better characterization of the instability. As pointed out 
in Paper I, where a mass ratio of 0.05 was also considered, 
mass ratios below 0.1 are more difficult to simulate, so we 
do not consider them in this study. According to the results 
of Paper I, they should lead to very similar results to the 
Md/Mbh =0.1 case. 
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Table 2. Initial models. The following parameters are listed: mass of the black hole MbHi disc-to-hole mass ratio Md/Mbhi slope of 
the specific angular momentum in the disc a, ratio of the potential barrier at the inner edge over the potential at the cusp AWin/| Wcuspl, 
mass flux in the stationary regime rhstat, minimum and maximum radii of the grid rmin and rmax, radius of the cusp Tcusp, radius of the 
centre rcontro (all these radii are in unit of the gravitational radius rg), and orbital period at the centre of the disc t^j.^ (in units of rg/c). 
The last column lists the timescale associated with the runaway instability as defined in Section l5.4l In all cases, the EoS parameters are 
K = 4.76 X 10^* cgs and 7 = 4/3. As mentioned in the text, models 3c' and 3c" differ only from model 3c by the value of the efficiency rj 
of the angular momentum transfer. 
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Notes: 

Initial value of the inner radius of the grid. When the spacetime is evolving, this radius increases with time (see text). 

(2) First value: inner part of the grid (high resolution); second value: outer part of the grid (low resolution). 

(3) This value is given for the case where both the mass and the spin of the black hole increase with time. 



5.3.2 Equation of state 



5.3.3 Mass flux 



We fix the adiabatic index to 7 = 4/3 and the polytropic 
constant to k = 1.2 x 10^^ Y^'^ with = 0.5. This cor- 
responds to an EoS dominated by the contribution of rela- 
tivistic degenerate electrons (the typical density in the disc 
is ~ 10^^-10^^ g/cm^). We note that such simplified EoS 
is nevertheless adequate f or our purpose since the work of 
iNishida fc Eriguchi lll99d) showed that, for stationary mod- 
els, the effects of realistic EoS on the stability of constant 
angular momentum discs is negligible, the discs being un- 
stable in all cases. 



The value of the potential barrier at the disc surface is 
strongly related to the acrretion mass flux in the stationary 
regime, as defined in Section |1] As we discussed in Paper I, 
this allows to limit the value of AWin in a domain where 
the corresponding mass flux is in the typical range expected 
for GRBs, 10~^-10^ Mqs~^, i.e. high values many orders 
of magnitude above the Eddington limit. Such very high 
mass fluxes are also found in the numerical simulations of 
compact binary coalescence or coUapsars mentioned above. 
In the models listed in Table |21 we have either fixed the 
potential barrier to a constant fraction of the potential at 
the cusp, AWin/l VKcuspl = 0.75 or 0.5, or we have adjusted 
it to have a constant mass flux in the stationary regime, 
1.3 Mqs"'^. The preparation of this last sequence of models 
is much more time consuming than for the rest of models 
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as it implies the following procedure: (i) construction of the 
initial model for a given first guess of AM^in ; (ii) evolu- 
tion of this model with the hydrodynamical code with fixed 
spacetime, up to the stage where the stationary regime is 
reached ; (iii) from the comparison of the mass flux mstat 
with the expected value, either start again at step (i) with 
a new guess of APFin or stop. 

All these considerations have led to the preparation of 
26 initial models listed in Table |21 They are splitted in four 
sets. Sets 'a' and 'b' correspond to a disc-to-hole mass ratio 
of 1, and sets 'c' and 'd' correspond to a mass ratio of 0.1. 
For each mass ratio, the two sets of models correspond to 
a different assumption for the potential barrier at the disc 
surface. The successive models within each set are obtained 
by increasing the slope a of the angular momentum distri- 
bution in the disc, starting from a = 0. 

5.4 Simulations 

Each of the 26 models listed in Table |5| is evolved three 
times: in the first series of runs (series one hereafter) both 
the mass and the spin of the initial (Schwarzschild) black 
hole are kept constant to their initial values. Correspond- 
ingly, in the second series of runs (series two hereafter) only 
the mass of the black hole is allowed to increase according to 
the procedure specified in the preceding section, while the 
black hole angular momentum is kept fixed to zero. Finally, 
in the third series of runs (series three hereafter) both the 
mass and the spin of the black hole are allowed to increase. 
The growth rate is monitored through the mass and angular 
momentum transfer from the disc to the black hole across 
the innermost radial zone. Keeping all these various possi- 
bilities into account gives a total of 78 simulations. These 
simulations are usually stopped at about t ~ 200 torb if the 
runaway instability has not occured before, which in two 
cases (models 3b and 4c) does not allow for the correct de- 
termination of the timescale of the instability (a lower value 
IS given m Table I^J, and in one case (model 3d) does not al- 
low to determine if the model is stable or not (a lower limit 
on the timescale of the instability in the latter case is also 
given in Table |5J. 

As mentioned in S ect ion f5 . 2 1 most of our initial models 
(of series three) are evolved assuming a conservative value 
for the efficiency of the angular momentum transfer, r; = 0.2 
(cf. Eq. H58|l ). In addition to these simulations we also evolve 
model 3c with two different values for this efficiency, namely 
»7 = 0.5 (model 3c') and rj — 1.0 (model 3c"). This allows us 
to check how sensitive our results on the runaway instability 
are on the amount of angular momentum transferred from 
the disc to the black hole. 

The time evolution of the mass flux for the difi'erent 
models is plotted in Figs. |2l and 1101 Additional time evo- 
lution plots for the mass of the disc and the mass of the 
black hole, as well as for the spin of the black hole, appear 
in Figs. 111! and 1121 respectively. For the sake of clarity in 
the presentation the time evolution plots of the mass flux in 
Fig. Elinclude all models listed in Table |21 (apart from mod- 
els 3c' and 3c"), but only for series three, i.e. when both the 
mass and the spin of the black hole increase with time. Cor- 
respondingly, the mass flux evolution displayed in Fig. 1101 
includes a representative selection of the models listed in 
Table |5| for all three possible combinations we consider for 



the evolution of the mass and spin of the black hole {dotted 
line: black hole of constant mass and angular momentum, 
series one; dashed line: black hole of increasing mass, series 
two; solid line: black hole of increasing mass and spin, series 
three). Similarly, the curves displayed in Figs. llll and ll2l cor- 
respond to all 26 initial models for evolutions of series three. 
In Figs.l^ nUllllI aud i 121 the left panels correspond to a disc- 
to-hole mass ratio Mu /Mbh = 1, i.e. set 'a' (top-left) and set 
'b' (bottom- left), and the right panels to Mo/Mbh ~ 0.1, 
i.e. set 'c' (top-right) and set 'd' (bottom-right). The evolu- 
tion of the mass flux in models 3c, 3c' and 3c" is plotted in 
an additional figure (see Fig. llBl belowl to illustrate the role 
of the efficiency rj of the angular momentum transfer. In all 
these figures the time is given in units of the orbital time at 
the centre of the torus, whose precise value can be found in 
Table |5| The results for each four classes of models listed in 
Table |5| are discussed in the following sections. 



5.4-1 Early time evolution of the mass flux 

The time evolution of the mass flux is presented in Figs. El 
and 1101 For all models and all series the early evolution 
is qualitatively identical: (i) initially the mass flux has a 
constant and very low value. This is due to the spherical 
accretion of the 'dust' solution which we have used to fill 
the region outside the initial torus. This component of the 
accretion rate is negligible in the subsequent evolution of 
the system; (ii) the mass flux increases abruptly and rapidly 
reaches a stationary regime (in about 1 — 2 orbital periods) . 
This is triggered by the fact that the initial torus overfills 
its Roche lobe, due to the positive potential barrier A Win. 
As shown in Paper 1 for the case of a Schwarzschild black 
hole and tested again in Section 01 for an a = 0.9 Kerr black 
hole, the value rristat of the mass flux in this regime follows 
closely rh oc AWt , which is the theoretical expectation for 
7 = 4/3 (Kozlows kTet alJll978l) . Tabled gives the value of 
the mass flux in the stationary regime for all models, which is 
independent of the considered series of runs. Notice that the 
potential barrier AVKin of the models of class 'd' have been 
adjusted to yield a constant mass flux rhstat = 1.3 Mqs"^ 
which is indicated in the bottom- right panel of Figs. |51 and 
IIUI by a horizontal dotted line. After this initial step, the 
time evolution of the mass flux differs strongly for the three 
series of runs. 



5.4.S The runaway instability of constant angular 
momentum discs 

As we showed in Papers 1 and 2, for constant angular mo- 
mentum discs (models la, lb, Ic, and Id; a = 0) the 
time evolution of the system changes dramatically when 
the spacetime dynamics is taken into account, i.e. when ei- 
ther Mbh only or Mbh and Jbh increase. In either case 
the runaway instability, reflected in the rapid growth of the 
mass accretion rate, appears on a dynamical timescale of 
~ 4.3 toib ~ 6.9 ms for model la (respectively, ~ 39 torb ~ 
88 ms for model lb, ~ 7.5 torb ~ 12 ms for model Ic, and 
~ 28 toib ~ 51 ms for model Id). The time derivative of 
the mass flux also increases, which implies the rapid diver- 
gence of 771. As already mentioned, this effect is more easily 
observed for models with Md/Mbh = 1. By comparing the 
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Figure 9. Time evolution of the mass flux: the mass flux is plotted as a function of time (normalized by the orbital time) for all 26 
models listed in Tabloid] in the case where both the mass and the spin of the black hole increase with time according to the procedure 
described in Section 15.21 For models Id to 6d (bottom-right panel), a horizontal dotted line indicates the constant value of the mass 
flux in the stationary regime. Unstable models are characterized by the sudden divergence of the mass flux in a few orbital periods. The 
transition between stable and unstable models is clearly visible in all four panels. 



case where only Mbh increases (dashed line in Fig. with 
the case where both Mbh and Jbh change (solid line), it 
can be seen that the instability appears slightly later in the 
latter case, i.e. the angular momentum transfer tends to dis- 
favour the insta bility. This is in agreement with the result 
first reported bv lWilsonI (^S^) for constant angular momen- 
tum tori, who showed using stationary models that when the 
spin of the black hole increases, the cusp moves toward the 
black hole and the torus is stabilized. 



5.4-3 Stabilizing effect of an increasing a 

The dynamical evolution of the disc, however, changes no- 
tably for non-constant angular momentum discs. The larger 
the slope a in the angular momentum distribution in the 
disc is, the more stable the models become, as it is visible 
in Figs.lUland llUI at least for the evolution times considered 
in our simulations (~ 200 forb)- This result applies to both, 
models of series two and series three. In the former case the 
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Figure 10. Time evolution of the mass flux: the mass flux is plotted as a function of time (normahzed by the orbital time) for a 
representative selection of the models listed in Tabloid The line styles indicate the three different cases considered: black hole of constant 
mass and angular momentum (dotted line), black hole of increasing mass (dashed line), and black hole of increasing mass and spin (solid 
line). For models Id, 2d, 4d and 6d, a horizontal dotted line indicates the constant value of the mass flux in the stationary regime. Notice 
that when the mass and the spin of the black hole are kept constant to their initial values the discs always remain stable irrespective of 
their angular momentum distribution. 



instability is supressed for larger values of a than in the lat- 
ter case. For instance it is possible to see in Fig. Cni that 
models la, 2a, 3a and 4a [a < 0.08) undergo the runaway 
instability in series two and three, model 5a {a — 0.085) 
is unstable only if the mass of the black hole increases (se- 
ries two, dashed line) and is stable when the increase of the 
black hole spin is also taken into account (series three, dot- 
ted line), and models 6a, 7a and 8a (a > 0.09) are stable in 



all series. The same transition appears in model 3b for set 
'b', model 5c for set 'c' and model 2d in set 'd'. 

Fig. |U] where all models of series three (both Mbh and 
Jbh increase) are plotted, shows clearly this transition be- 
tween unstable and stable models. Correspondingly, Fig- 
ure 1131 shows the location of these models in the mstat -a 
parameter plane. Unstable models are indicated with a black 
circle and stable models with an empty circle. The region of 
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Figure 11. Time evolution of the mass: the mass of the disc (solid lines) and the mass of the black hole (dotted lines) are plotted 
as a function of time for all models listed in Table \^ only in the case where both the mass and the spin of the black hole increase with 
time. For models with a disc-to-hole mass ratio A/d/A/bh = 0.1 (right panels) the corresponding insets show a blow up of the evolution 
of the mass of the black hole in the latest stages of the simulation. 



stability appears clearly in the top-left corner of the figure. 
The critical slope cvcr separating the unstable and stable 
regimes depends on the initial mass flux. It is plotted in 
Fig. I13l and well fitted (but only for four sets of models!) by 
the relation Ocr ~ 0.018 -I- 0.021 (mstat/M©). The frontier 
between the two regions in the riistat-a plane is plotted us- 
ing this simple fit. Notice that this critical slope is always 
well below the Keplerian limit a = 0.5, and even well below 
the typical slope found in numerical simulations of compact 
binary coalescence and coUapsars. For very high initial mass 



fluxes (rhstat 3> 10 Mqs"^), the extrapolation of our esti- 
mate of the critical slope acr indicates that models become 
more unstable. However, the discussion of the runaway in- 
stability is probably not relevant for such high mass fluxes as 
the stationary value of the mass flux is already large enough 
to destroy the disc in a few orbital periods. 

The low value we find for the critical slope Qcr 
(ocr G [0.085; 0.09] for the set of models 'a'; respec- 
tively, [0.015; 0.02] for set 'b', [0.055; 0.06] for set 'c' and 
[0.04; 0.055] for set 'd') indicates that a small increase of 
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Figure 12. Time evolution of the spin of the black hole: the Kerr parameter a is plotted as a function of time for all models listed 
in Table 121 Unstable models with disc-to-hole mass ratio A/d/A/bh = 1 (left panels) spin up the initial Schwarzshild black hole up to 
o/A/bh ~ 0.3 — 0.35. Correspondingly, unstable models with disc-to-hole mass ratio A/d/A/bh = 0.1 only reach a/A/BH ^ 0-04 — 0.05. 



the specific angular momentum with the radial distance 
strongly stabilizes the disc, completely suppressing the run- 
away instability. This stabilizing effect that we find in our 
time-dependent general relativistic simulations is in com- 
plete agreement with previous studies based on station- 
ary sequences of equilibrium configurations, either using a 
pseudo-Newtonian potential (DaiEnc & Mochkovitch 19971) 
or in general relativity llAbramowicz et al.ill998i) . Moreover, 
the precis e value of the critical slope Qcr is also in good ac- 
cordanc e. iDaigne fc M ochkovitcbl lll997l) (Newtonian) and 
lAbramowkz^^dTi^g^ (relativistic) found a^r ~ 0.07 and 



Ocr = 0.14, respectively, for a = 0, Mbh = 2.44 Mq and 
Md/Mbh = 0.148 (k and 7 being the same than the ones 
we use in this paper). The dependence of acr on the mass 
fiux is of course masked in these previous calculations as 
they assume an instantaneous transfer of mass. 

5.4.4 Time evolution of the disc and the black hole mass 

In Fig. llll we plot the time evolution of the mass of the disc 
(solid lines) and the mass of the black hole (dotted lines) for 
all models in series three (both Mbh and Jbh increase). The 
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Figure 13. Left: the best estimate of tlie critical slope ctcr is plotted as a function of the initial mass flux rhstat- The four points 
corresponds to the mean value between models 5a and 6a (set 'a'), models 3b and 4b (set 'b'), models 4c and 5c (set 'c'), and models 
2d and 4d (set 'd'). The dashed line indicates the empirical relation acr — 0.018 + 0.021(m/MQS~^). Right: all models (series three) are 
plotted in the rristat-Cf plane. Filled circles indicate unstable discs and empty circles stable discs. The dashed line separating the stable 
and unstable discs corresponds to the empirical relation shown in the left panel. 



unstable behaviour of models la to 5a, lb to 3b, Ic to 3c, 
and Id to 3d is characterized by the sudden decrease of Md 
and the corresponding increase of Mbh. In the right panels, 
in which the discs are much less massive than in the left 
panels, it is seen that the high efficiency of the disc-to-hole 
mass transfer does not imply an important growth of the 
black hole mass. The inset of these two panels focus on the 
late time evolution when the increase is more pronounced. 

In all unstable models, at least 93 % of the disc mass has 
been accreted by the end of the simulation (not displayed 
in the figure). The minimum value is reached for model 2d 
(Mbh = 2.5 M© and A/d = 0.25 M©), in which the final 
mass of the disc is only 0.016 Mq. High initial mass ratios 
lead to more extreme situations: for instance in model la 
(Mbh = Md = 2.5 Mq), the final mass of the disc is only 
2.8 X 10~* M0 and the final mass of the black hole is 5 M©. 
On the other hand, for models 6a, 7a and 8a {a — 0.09, 0.1 
and 0.15), 4b, 5b and 6b (a = 0.02, 0.025 and 0.05), 5c, 
6c and 7c (q = 0.06, 0.07 and 0.075), and 4d and 5d (a = 
0.055 and 0.06), the stability properties of the non-constant 
distribution of angular momentum are reflected in the small 
loss of the mass of the disc throughout the evolution: less 
than 3.7 % of the disc mass has been accreted in sets 'a' 
and 'b' (maximum reached in model 6a) and less than 27 
% in sets 'c' and 'd' (maximum reached in model 5c). The 
higher fraction found for the case of a low disc-to-hole mass 
ratio Mo/Mbh ~ 0.1 is explained below in the subsection 
devoted to the long-term evolution of stable discs. 

As in Paper I we define the timescale of the runaway 
instability tmn as the time needed for half of the disc mass 
to fall into the hole. An interesting outcome of the present 
study is to provide the first numerical estimation of this 
timescale for non-constant angular momentum discs. The 



last column of Table |21 hsts tmn (in units of the orbital 
period) for all models (series three). It can be as short as 
trun ~ 4.3 torh ~ 6.9 ms (model la) and is always shorter 
than 1 s in all models considered in our sample. The max- 
imum values are the lower limits obtained for model 4c, 
trun > 200 toih ~ 0.38 s, and model 3d, Uun > 600 torb ~ 
0.6 s, the latter case being possibly stable. 

Fig. ll4| plots the dependence of the timescale of the run- 
away instability with a (left panel) and with rhstat (right 
panel) - which is equivalent to AVKin as explained above - 
for the four sets of models in series three. In order to identify 
more clearly the effect of each one of these two parameters 
we have isolated in Fig. ll5l those models that allow to study 
the dependence of the timescale with a for constant rhstat 
(left panel) and with rhstat for constant a (right panel). As 
for the case of constant angular discs (see Paper I), we find 
now that for non-constant angular momentum discs the run- 
away instability occurs also faster when the initial mass fiux 
in the stationary regime is larger. As plotted in Fig. ll5l (right 
panel), we get an empirical dependence tmn oc m'^^lf for 
a — and an even steeper decay tmn oc rn~^l^f for a = 0.05. 
Notice that the slope -0.5 found for q = is slightly smaller 
than the value -0.9 found in Paper I. This is due to the 
stabilizing eflect of the angular momentum transfer which 
was not taken into account in Paper I. On the other hand, 
the timescale of the runaway instability increases exponen- 
tially with a until it is finally suppressed when the critical 
slope Qcr is reached. The left panel of Fig. 1151 shows that, 
in the case where riistat ~ 1.3 Mq, the timescale follows the 
empirical dependence tmn ~ 50 exp (a/0.021) ms. This ex- 
ponential increase shows again how strong is the stabilizing 
effect of a non-constant angular momentum distribution in 
the disc. 
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Figure 14. Timescale of the runaway instability: all models. Left: the timescale iiun is plotted as a function of the slope of the 
angular momentum distribution. Right: the same timescale is plotted as a function of the initial mass flux rhstat which is induced by the 
potential barrier AWin at the surface of the disc. 
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Figure 15. Timescale of the runaway instability: selected models. Compared to Fig. 1141 models with a constant mass flux have 
been selected and included in the left panel, and models with constant angular momentum slope have been selected for the right panel. 
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On the other hand Fig.|n|shows (compare set 'a' top-left 
and set 'b' bottom-left) that a smaller overflow of the initial 
marginally stable potential barrier (resulting in a smaller ini- 
tial mass flux) tends, not surprisingly, to defer the occurence 
of the instability. The timescale tmn is typically ten times 
longer in set 'b' compared to set 'a'. For the same reason 



Fig.|5|aIso indicates that the smaller the initial mass flux is, 
the smaller the slope a necessary to suppress the instability 
is. The critical slope is ~ 0.09 for set of models 'a' and only 
0.02 for set 'b'. Finally, Fig. [^illustrates that in comparison 
with the strong dependence on rhstat and a, the dependence 
of the runaway instability on the disc-to-hole mass ratio is 
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rather weak (compare for instance set 'a' (top-left panel) 
and set 'c' (top-right panel) in Fig. (nj. 



5.4.5 Time evolution of the black hole spin 

The transfer of mass and angular momentum from the disc 
to the black hole leads to the gradual increase of the rotation 
law of the initially non-rotating black hole. Fig. ll2l shows the 
time evolution of the black hole spin q/Mbh for all models 
considered. Once again there is a clear distinction between 
stable and unstable situations. The stable discs only transfer 
a very small fraction of their angular momentum to the black 
hole. At the end of the simulation the black hole spin is less 
than a/Afan = 1-5 x 10~^ (maximum reached for model 5c) 
and can be as low as 1.9 x 10"'' (model 6b). On the contrary, 
for the unstable discs the transfer of angular momentum is 
considerably higher, accelerating rapidly during the runaway 
instability. As a result, in less than 100 orbital periods the 
initial Schwarzschild black hole turns into a mildly rotating 
Kerr black hole with q/Mbh ~ 0.3 — 0.35 for sets 'a' and 'b' 
(Md/-A/bh = 1) and into a slowly rotating black hole with 
a/MBH<0.05 for sets 'c' and 'd' (Md/Mbh = 0.1). In the 
second case, the small increase of the spin is directly related 
to the fact that for such disc-to-hole mass ratio the reservoir 
of angular momentum is much smaller. 

These results on the spin of the black hole have been 
obtained for a conservative value of the efficiency 77 = 0.2 
of the disc-to-hole angular momentum transfer. Figure 1161 
illustrates the dependence of our results with -q. The mass 
flux is plotted as a function of time for models 3c, 3c' and 
3c". This figure clearly shows that increasing the efficiency 
of angular momentum transfer has a stabilizing effect. While 
model 3c is unstable for r) — 0.2, model 3c' is stable for 77 = 
0.5 and model 3c" is even more stable for rj — 1. Therefore, 
except for very low values of the slope a, a high efficiency 
?7 would probably not lead to more rapidly rotating black 
holes but most likely to stable discs. On the other hand, 
if the initial black hole is already rotating, as it is inferred 
from the existing numerical simulations of compact binary 
coalescence or coUapsars, this effect might result in black 
holes with high values of a/MBu. This will be studied in a 
forthcoming paper. 



5.4.6 Long-term evolution of stable discs 

The long-term behaviour of the mass flux evolution in stable 
discs simply reflects the fact that the initial accretion rate 
is far from zero but rather highly super-Eddington. For in- 
stance in the set of models 'a', the mass flux is ~ 3.2 Mq/s 
for model 6a, ~ 2.5 Mq/s for model 7a, and ~ 0.69 Mq/s 
for model 8a. Hence, by integrating m along the entire time 
of the simulation, it is possible to check that the total mass 
transferred to the black hole becomes, asymptotically, of the 
order of the initial mass of the disc, which results in the end 
of the accretion process. This effect is more pronounced for 
smaller disc-to-hole mass ratios (sets 'c' and 'd') as the accre- 
tion timescale ~ rhatnt/Mo is shorter, or in other words, as 
the reservoir of mass is smaller. Furthermore, as we already 
noticed in Paper 2, a second effect is the gradual decrease of 
the mass flux during the long-term evolution of stable mod- 
els (and of unstable models before the instability sets in). 




0.01 0.1 1 10 

Time t / t„. 



10= 



Figure 16. Time evolution of the mass flux: the mass flux 
is plotted as a function of time (normalized by the orbital time) 
for models 3c, 3c' and 3c", which correspond to three different 
values of the efficiency r] of the disc-to-hole angular momentum 
transfer: t] = 0.2 (conservative value used in all other simulations, 
solid line), t] = 0.5 (dashed line) and rj = 1 (dotted line). Notice 
the stabilizing effect of an increasing rj. 



This is simply due to the fact that for non-constant angular 
momentum discs the accretion process makes the cusp move 
faster towards the black hole than the inner radius of the 
disc (which is precisely why the instability is physically sup- 
pressed). Therefore, the potential barrier at the inner edge 
AVKin decreases with time and, consequently, also the mass 
flux. For this reason, as shown in Fig. |^ accretion becomes 
increasingly slower, which, in turn, makes numerically dif- 
ficult to follow the evolution until the entire disc has been 
fully accreted. 



6 CONCLUSIONS 

We have presented a comprehensive numerical study 
aimed at clarifying the likelihood of the so-called run- 
away instability in axisymmetric thick discs orbiting around 
(initially non-rotating) black holes. Our main conclu- 
sion, already pointed out in a preliminary investigation 
jFont fc D aign(j"2002cl) is that the runaway instability of 
(non self-gravita ti ng) thick d iscs around (rotating) black 
holes llAbramowicz et al.lll983^ is strongly suppressed when 
the distribution of angular momentum in the disc is assumed 
to be non constant. We have reached this conclusion using 
a numerical framework based on time-dependent hydrody- 
namical simulations in general relativity. This approach de- 
parts from the most common procedures followed in previ- 
ous studies of the runaway instability in the literature, which 
are generally based on the a nalysis of stationary models (see. 
e.g. lFont fc Daignd ^l2002b^ and references there in) . Despite 
the novel approach it is remarkable that our main result is 
in complete agreement with such previous studies based on 
stationary sequences of equilibrium c onfi gurations, either us- 
ing a pseudo-Newtonian pote ntial jPai gne & MochkovitchI 
Il997tl or in general relativity (lAbramowicz et al.lll998f) . We 
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further notice that our finding is also in good agreement with 
recent pertur bative analysis of ve rtically integrated discs 
performed bv lR.ezzolla et aO l)2003h . 

The present paper extends in a number of ways the 
preliminary results rep orted in our previous investigation 
llFont fc Daignd jiooi^). In the first half of the article 
(cf. Section we have shown how to build a series of 
thick accretion discs in hydrostatic equilibrium around a 
Schwarzschild or Kerr black hole. All relevant aspects of 
the theory of stationary thick discs around rotating black 
holes have been presented in great detail. In particular our 
description has been done for discs in which the specific 
angular momentum (per unit mass) is assumed to increase 
outwards with the radial distance according to a power law 
I — Kr°' . All possible equipotential configurations for a Kerr 
black hole (including the limiting case of a Schwarzschild 
black hole) have been discussed in detail for both, pro- 
grade {K > 0) and retrograde {K < 0) rotation, and sub- 
Keplerian (0 < a < 1/2), Keplerian {a — 1/2), and super- 
Keplerian (1/2 < a < 1) discs. The various possible geome- 
tries of these configurations have been summarized in Ta- 
bleland in Figs.lHtoEl These figures and table extend ear- 
lier results presented in lFont &c Daign c (2002b) where only 
constant angular momentum discs around Schwarzschild 
black holes were considered. 

The second half of the article has been devoted to pre- 
senting the results of a comprehensive set of models that we 
have evolved with our numerical code. The main result of 
these simulations has been to show, quite generically, that 
a small increase outwards of the specific angular momen- 
tum distribution in the discs, well below the Keplerian limit, 
results in a strong stabilizing effect. In particular we have 
demonstrated that the case a = (constant angular momen- 
tum disc) is always unstable when the black hole dynamics 
is taken into account, for both, Md/A/bh = 1 and 0.1, ir- 
respective of the value of the potential barrier at the cusp 
(slightly overfiowing the corresponding Roche lobe) and of 
the initial mass flux in the stationary regime. We note in 
pa ssing that a sim ilar conclusion has been recently obtained 
bv lZanotti et al ] 120031 for a Schwarzschild black hole and 
for different initial data, in which models that initially are 
in strict stable equilibrium are conveniently perturbed to 
trigger the accretion process. However, we have shown that 
discs with small non-zero values of a are stable, and we 
have provided an estimate of the critical value of a sepa- 
rating stable and unstable models, for various disc-to-hole 
mass ratios and mass fluxes in the stationary regime. More 
precisely we have found that for Mb/Mbh = 1, the critical 
slope is Ucr ~ 0.085 — 0.09 for rhstat ~ 3.2 — 3.6 Mqs~^ and 
Qcr ~ 0.015 — 0.02 for mstat ~ 0.1 Mqs"^. Correspondingly, 
for AIo/Mbh = 0.1, the critical slope is Qcr ~ 0.055 — 0.06 
for rhstat ~ 1.7 — 1.9 Mqs"^ and Qcr ~ 0.05 for riistat ~ 
1.3 Mqs^. These values are slightly smaller than the values 
found in stationary studies (0.14 in the relativistic case for a 
disc-to-hole mass ratio of ~ 0.15 jAbramowicz et 'al]ll998^ 'l. 
However these studies assume an instantaneous mass and 
angular momentum transfer from the disc to the black hole, 
and our study clearly indicates that the critical slope ftcr in- 
creases with the mass flux rhstat . For the first time we have 
shown the suppression of the runaway instability in non- 
constant angular momentum discs with a time-dependent 
relativistic calculation, and we have been able to estimate 



the lifetime of thick discs orbiting a black hole in both, the 
unstable and stable cases. The very strong stabilizing effect 
of the non-constant angular momentum distribution in the 
disc is demonstrated in our simulations by the exponential 
increase of the timescale for the instability to set in with the 
slope a. 

The black hole plus thick discs systems we have sim- 
ulated with Md/A/bh = 0.1 are very close to what 
numerical simulati ons of compact bina r y mergers and 
coUapsars predict jRuffert fc Jantol |l999l: I Kluzniak fc Led 
| l998t iMacFadven fc WooslevI 119981 : IShibata fc Urvull200d : 
[^ove^jTl200(ll . Therefore, the results regarding this sub- 
set of models should be considered as, perhaps, much closer 
to what may be happening in realistic discs. However, even 
in this case, three important simplifications have been made: 
(i) the black hole is initially non-rotating whereas both in 
coUapsars and binary mergers it is expected to form with 
an initial spin parameter a/AfBH^O.5; (ii) the self-gravity 
of the disc is not included; and (iii) magnetic fields, which 
are believed to play an important role on the dynamics of 
those objects, are also not included. 

Certainly, along with analyzing the effects of varying 
the EoS of the matter in the disc, these are all issues that 
we plan to consider in the near future in order to extend 
the present study. Among those, the inclusion of an appro- 
priate treatment to account for the self-gravity of the disc 
is the most dif ficult task despite its obvious r elevance. Pre- 
vious studies (iKhan na. fc Chakrabartil Il992l : iNishida et all 
119961: 1 Masuda 6^01^9^ . based on either time-dependent 
simulations with pseudo-Newtonian potentials or station- 
ary models in general relativity, indicate that self-gravity 
seems to have a strong d estabilizing effect. It has been ar- 
gued jNishida et al.lll996^ that in self-gravitating accreting 
tori the cusp is pushed towards the black hole by the grav- 
itational force of the discs which, in principle, should act 
against the instability. However, during the accretion pro- 
cess, the interplay between the gravitational forces of the 
black hole (increasing) and of the disc (decreasing) make 
the cusp move closer to the centre of the torus than in 
non-self-gravitating models, thus favouring the instability. 
This important result is yet to be proved by means of time- 
dependent simulations in full general relativity. Unfortu- 
nately, successful attempts to long-term (numerically) sta- 
ble simulations of thick discs orbiting black holes by solving 
the coupled system of Einstein and general relativistic hy- 
drodynamic equations are still out of reach of p resent day 
numerical relativity codes (see, e.g. iFond i2003l) and refer- 
ences there in for an up-to-date list of such codes). Such 
simulations should provide a definite answer to the issue of 
whether the runaway instability exists or not. 

We notice that neglecting the self-gravity of the disc is a 
more valid assumption the smaller the disc-to-hole mass ra- 
tio becomes. Therefore, with respect to the models analyzed 
in this paper we expect that for discs with Mu/Mbh ^0.1 
the strong stabilizing effect of non-constant angular momen- 
tum distribution demonstrated in our study could still over- 
come the destabilizing effect of the much less relevant disc's 
self-gravity. 

The incorporation of magnetic fields seems to be 
more easil y conceivable. Recen tly, iDe Villiers fc HawlevI 
(2003) and ICammie et all (|2QQ^ have presented magneto- 
hydrodynamical simulations of constant angular momentum 
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tori in general relativity. These simulations assume much 
smaller accretion mass fluxes than in the present study and 
do not include any time evolution of the spacetime. However, 
they indicate a promising direction to include the effect of 
magnetic fields in future simulations of the runaway insta- 
bility and hence suppress the third limitation listed above. 

Nevertheless, the first aspect that we plan to consider in 
a forthcoming study is the time evolution of discs orbiting 
black holes which are initially rotating. This will allow us 
to pinpoint the dependence of Ocr with a. We notice that 
our hydrodynamics code has already proved in this study 
its ability to deal with the Kerr spacetime, and the method 
of construction of the initial state in the case of a rotating 
black hole has been fully described in this paper. Therefore, 
our upcoming study does not involve any new numerical 
development. 



APPENDIX A: COEFFICIENTS OF THE KERR 
METRIC AND THEIR DERIVATIVES IN THE 
EQUATORIAL PLANE 

In this appendix we provide the expressions for the Kerr 
metric coefficients and their derivatives in the equatorial 
plane, which are needed in the construction of the equilib- 
rium tori. In the following we assume that the black hole 
mass is Mbh = 1 (therefore < a < 1) and we use the 
standard Boyer-Lindquist (spherical) coordinates {t, r, 8, (j)) 
to describe the Kerr metric (cf. Eq. Q in Section I^J: 

2r 

9tt = "1 + -r— T-:7777T (Ai) 



9t:i> = 



+ cos^ 9 

2ar siv? 6 

+ cos^ 9 

2 2a'^rs\v? 9 2 \ . 2, 

■r + — 5 5-5 + a I sm ■ 

-I- a-^ COS'' 9 

-f cos^ 9 



— 2r + 

2,2 2 r 
r -f a cos 6 
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(A3) 

(A4) 
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The pseudo-distance vj, the radius of the horizon rh and the 
radius of the ergosphere are given by: 



Txj'^ = gi^ - gttg4,<t> = (r^ - 2r + a^) sin^ ( 



,{9) = 1 + 



(A6) 
(A7) 
(A8) 



The metric coefficients in the equatorial plane which are 
needed in Section |5| are then given by: 
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The first and second derivatives of these coefficients are also 
needed: 
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The pseudo-distance to the axis, Eq. IIA6L is also a useful 
quantity. We give here its expression as well as its derivatives 
in the equatorial plane: 



vj"^ = -2r + 
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dr 
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APPENDIX B: VON ZEIPEL'S CYLINDERS 
Bl Equation of the cylinders 

In this appendix we describe the structure of the von Zeipel's 
cylinders (i.e. the surfaces of constant angular momentum 
per unit inertial mass and constant angular velocity, which 
coincide for a barotropic equation of state) in the Kerr space- 
time. We use the same notations and conventions than in 
Section |5| and Appendix El We limit the discussion to the 
region outside the black hole horizon and to those cylinders 
which have an intersection with the equatorial plane, where 
we assume that the distribution of angular momentum per 
unit inertial mass lcq{r) is given. The equation of the von 
Zeipel cylinder corresponding to I = Iq is: 

gttlo + gt,t> {1 + ^olo) + g<i,4,^o = , (Bl) 

where fig is the angular velocity corresponding to Iq. Both 
quantities are related by 

9t4> + 9ttl 



g,i><i> + gt<t>i 



(B2) 
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This leads to the equation of the von Zeipel's cylinder inter- 
secting the equatorial plane at radius ro: 

= Z{r, 9 ■,ro) = [gttgt<t>{ro) ~ gt,pgtt{ro)] llq{ro) 
+ [gttg<i,4.{ro) - g<l><t>gtt{ro)] icq(r-o) 
+ [gt4,g4><t>{ro)- g4,4>gt4,{ro)] . (B3) 



B2 Intersection with the equatorial plane 

We investigate now the structure of the von Zeipel's cylin- 
ders in the vicinity of the equatorial plane. For 6 = 7r/2, 
2(ro,7r/2 ;ro) = and the derivatives of Z{r,9 ;ro) at 
(ro, 7r/2) can be computed from the coefficients of the metric 
and their derivatives (see Appendix IXt: 



— (ro,7r/2 ;ro) 

— (ro,7r/2 - ro) 
— {ro,n/2 -ro) 
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(B4) 
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drde 

O'^ Z 2 

-^(ro,7r/2 ;ro) = (ro - 2ro + a'^) Q (ro, «cq(ro)) , 

(B8) 

where the two polynomials V and Q are given by 

P(ro,/o) = all + (rl - 3rl - 2a^) lo + a{3rl + a^) {B9) 

Q{ro,lo) = 2all- {rl+4a'')lo + 2a\ (BIO) 

Then, in the vicinity of the equatorial plane, for r = ro + dr 
and 9 = n /2 + d9, the equation of the cylinder is given by 

= Z{ro + dr,-K/2 + d9 ■,ro) 
2 

= — {ro,laqiro)) dr + 

\ Q (ro, /cq(ro)) [dr^ + (ro - 2ro + a) d9^] + ... 

(Bll) 

There are two possible cases: 

(i) In the normal case, V {ro,lcq{ro)) 7^ and Eq. ijBllfl 
implies that the first derivative dr/d9 along the cylinder in 
the equatorial plane vanishes. In the vicinity of (ro,7r/2), 
the equation of the cylinder becomes 



r = ro 



ro - 2ro + Q{ro,kci{ro)) 



2ro 



(B12) 



It means that the von Zeipel's cylinder starts tangential to 
the circle of radius ro and moves aside from this circle with 
an increasing (respectively, decreasing) radius for Q/V < 
(respectively Q/V > 0) at (ro, Zcq(ro)). 

(ii) In the critical case where P (ro, /eq(»'o)) = 0, 
Eq. I|B11^ implies that the cylinder is divided in two 
branches which intersect in the equatorial plane. The inter- 
section (ro,'7r/2) appears as a cusp. The equation of these 
two branches in the vicinity of the cusp is simply given by 



: ro ± V^o - 2ro + a' (^9 - |) 



(B13) 



In the case where a > (rotating black hole), the two 
polynomials V and Q can be studied as a function of lo for a 
fixed radius ro. Outside the horizon V has two positive roots 
for rh < ro < r^^^ and two negative roots for ro > rp^'. It is 
negative if lo is between the two roots when they exist and 
positive everywhere else. The radii and are the roots 



of rg - 6rg + 9ro - 4o^ with r^ < r^^'' < 3 and 3 < r^"^ < 4 



For a — 0, r-g^^ — r\^' = 3 and for a = 1, r^''' = r^ = 1 
and r'^'^ = 4. The radius r^^^ crosses the ergosphere for 
a = \/2/2. The behaviour of Q is much simpler: it has always 
two positive roots and it is negative when Iq is between these 
two roots. 

The roots of V (respectively, Q) are plotted as a thick 
(respectively, thin) line in the ro-Zo plane in Fig. IBll The 
region where Q/V < appears shaded in the plot (criss- 
crossing lines). In the same figure, the von Zeipel's cylinders 
are plotted for a constant angular momentum /cq(?'o) = '0 
with different values of |/o|- The corresponding values of the 
angular momentum are indicated as an horizontal thick line 
in the ro-/o plane in each plot. In each case, the x-y plane 
has been divided in two regions: for y > 0, the cylinders 
are plotted for lo = +|^o| (prograde orbits) and for y < 
for lo = —\lo\ (retrograde orbits). When the line /oq = 'o 
intersects the curve V{ro,lo) ~ 0, it is indicated by a thick 
dot in the ro-/o plane and the corresponding two branches 
critical cylinder is plotted in the x-y diagram with a thick 
line. Notice that for a = 1 the only possible critical cylinder 
corresponding to prograde orbits is obtained with ro = rh = 
1 and Zo = 2. Its equation is given by (r — 1)^ cos'* 9 — (r** — 
4r -I- 3) cos'^ 9 — (r* — 3r'^ -I- 2r) = (last column, third row) . 

In the particular case a = (non-rotating black hole), 
the equation of the von Zeipel's cylinders becomes indepen- 
dent of the distribution of the angular momentum per unit 
inertial mass: 



= (ro - 2)r^ si: 



2)r^. 



(B14) 



There is always one critical cylinder with a cusp located at 
ro = 3. This case is also plotted in Fig. IBll (first column). 



APPENDIX C: DERIVATIVE OF THE 
POTENTIAL IN THE EQUATORIAL PLANE 

CI Mathematical expression in Kerr spacetime 

The derivative uicq(r) — dWoq/dr given in Eq. 1321 can be 
expressed from Eq. ijAB^ and the expression of the metric co- 
efficients and their derivatives in the equatorial plane given 
in Appendix 1X1 as 



1 All^ + 2BZeq + C 



2u}^ " gttHq + 2gt4,loq + g^,^, 



(CI) 



where the coefficients A, B and C are given by 

dgt4> , -2 dg^^ 
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2(r^ - 4r^ -f 4r - a^) 
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Figure Bl. The structure of the von Zeipel's cylinders outside the horizon: each column corresponds to a different value of 
the spin of the black hole a. The first row shows in the ro-lo plane the curve V = (thick line) and the curve Q = (thin line) as well 
as the region where Q/V < (shaded area), which allows to predict the behaviour of any cylinder in the vicinity of the equatorial plane 
(see text). The interior of the horizon has been shaded with crossed lines. For each value of a, the von Zeipel's cylinders arc shown in 
the x-y plane for ioq(''o) = iO.5 for a = 1/3 (resp. ±1 for a = 2/3 and a = 1) in second row, for icq(''o) = i2 in the third row and for 
leq{ro) = ±16 for a = 1/3 (resp. ±8 for a = 2/3 and a = 1) in the fourth row. The corresponding distribution leq{ro) are indicated in 
the first row as an horizontal thick line. In the case where there is a critical cylinder, it has been plotted with a thick line is the x-y plane 
and the corresponding (ro,/o) values have been indicated with a big dot in the rg-lo plane. For the second, third and fourth rows, the 
x-j/ plane is divided in two regions ■ y > corresponds to prograde rotation (ieq > 0) and y <0 to retrograde rotation (ieq(7'o) < 0). In 
the first column (a = 0, non-rotating black hole), the cylinders are independent of the distribution ieq(''o)- 



C - - [g^^ — - 2gt4.g4.4, — + gt^ j 

2(r^ + 2aV-4aV + a^) 
^2 • 



C2 Condition to have a finite value 

In the denominator of Weq(r), the first quantity 2w^ is ev- 
erywhere positive for r > rh and vanishes for r = rn, so 
that the derivative of the potential WccSj-) is infiuitc at the 
horizon. The second term is a polynomial of degree 2 for the 
variable Zeq- Its discriminant is simply tu^ so that there are 
always two distinct real roots outside the horizon: 
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Figure CI. The behaviour of 'ii!eq(r') = dWcq/dr outside the horizon: in the plane i — Isqir), the critical value l-j^{r) (respectively, 
lcr^{r)) where the derivative of the potential Weq{r) vanishes (respectively, becomes infinite) is plotted in solid line (respectively, dashed 
line). The interior of the horizon is the region shaded with crossed lines. The two radii rc = 2 and are indicated by a vertical dotted 
line. The area shaded with oblique lines corresponds to the region where the derivative Wcq(r) is negative. The left, middle and right 
panels correspond, respectively, to three representative values of the spin of the black hole: a = (non rotating black hole), a = \/5/3 
and a = 1 (maximally rotating black hole). Notice that the positions of the centre and the cusp in the disc are given by the intersection 
of the assumed rotation law /oq{»") with l-^{r), if any. 



gtt 



(C2) 



Then, Waqir) has always a finite value in the equatorial plane 
outside the horizon except for Zoq = itii^)- Notice that at 
the horizon lti{rh) = ±oo for a = and lti{rh) = Sr^/a > 
otherwise. Concerning the sign of the denominator of Wcq{r), 
there are two cases: (i) inside the ergosphere (r < 2) we 
have lci{r) > itii^) > and the denominator is negative for 
Icrir) < lcq{r) < Icrir)', (u) outside the ergosphere (r > 2), 
we have Z^(r) < < Z^(r) and the denominator is negative 
for lcq{r) < lcv{r) or lcq{r) > /i-(r). 

In Fig. Ell we plot l^{r) and /±(r) for a = (left), 
a = VE/S (middle) and a = 1 (right), indicating the region 
where Wcq{r) is negative. 



C3 Condition to vanish 

The discriminant A — — AC of the numerator equals 

^ _ 4{r^ ~2r + a^f 
r 

and is always positive so that w^^ij) can vanish under the 
condition /oq = {—B ± \/A)/j4. This leads to 

-a(3r^ -4r + a^) ± (r^ - 2r + a'^)r^ 



Zoq(r) = 



= ± 



(r^ + a^) =F 2ayF 



(r-2)Vr±a 

From equation I17L this condition is nothing else than 

^cq(r)=;^(r). (C3) 

Notice that at the horizon ij^(r'h) = ±oo for a = and 
^K(''h) = 2rh/a = ltv{rh) > otherwise. The coefficient A 



can be written as 
2 

r 



A 



((r - 2)v/f - a) ((r - 2)V? + a) 



(C4) 



so that outside the horizon, A is positive for r < r~ and 
negative otherwise. The radius r~ > 2 is the unique root 
for r > r-h of (r — '2)^/r — a. Then there are two cases 
for the sign of the numerator of u;oq(''): (i) for r < r~, we 
have l^{r) > l^{r) > and the numerator is negative for 
l^{r) < leq{r) < l^{r); (u) for r > r~, we have l^{r) < < 
l^{r) and the numerator is negative for lcq{r) < l^{r) or 
Zoq(r) > l+{r). 
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